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primarily  the  work  of  Keith  McLaughlin  and  forms  part  of  his  PhD  dissertation. 

Section  II  describes  the  analysis  of  array  data  recorded  1.9  km  from  the 
explosion  Liptauer  in  Yucca  Valley  of  the  Nevada  Test  Site.v  Coherency  declines 
gradually  with  both  inter-station  spacing  and  frequency.  Broadband  correlation 
across  the  400  m  array  is  greater  than  78%  for  the  horigotrlfal  components  and 
greater  than  64%  for  the  vertical  components^ 

■ >Section  III  presents  a  new  way  of  looking  at  frequency-wavenumber  spectral 
estimation  with  array  data.  .  The  differences  between  the  conventional 
beamforming  method  of  estimation  and  high-resolution  method  of  estimation 
is  quite  clear  when  ejcgsoa^ed  in  terms  of  the  eigenvalues  of  the  cross 
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^Section  IV  examines  the  problem  that  exists  when  the  velocity  structure 
is  complicated  to  the  extent  that  it  can  be  considered  to  have  a  random 
component.  .  This  gives  rise  to  a  random  component  on  the  seismogram  which  is 
usually  referred  to  as  signal-generated  noise  or  coda.  This  random  component 
can  be  separated  from  the  deterministic  part  of  the  seismogram  and  treated 
as  the  convolution  of  the  source  with  a  stochastic  Green's  function.  Estimates] 
of  this  stochastic  Green's  function  and  its  variance  are  useful  in  inferring 
bias  and  uncertainty  in  estimated  source  functions. 

'^Section  V  treats  the  problem  of  scattering  of  elastic  waves  by  small 
inhomogeneities.  The  solution  is  expressed  in  terms  of  a  moment  tensor 
expansion  of  the  properties  of  the  scatterer.  This  approach  is  convenient 
for  examining  the  trade-offs  between  shape,  fJeierogeneity ,  and  anisotropy  of 
scatterers.  It  reveals  that  a  general  scattere\can  not  be  modeled  by  a 
homogeneous  scatterer.  \ 

Progress  on  element  5  of  the  research  program  is  described  in  section  VI. 
In  order  to  provide  better  access  to  archived  seismic  data  and  facilitate 
computations  with  these  data,  a  Computational  Center  for  Seismology  (CCS) 
has  been  established.  The  organization  and  initial  efforts  of  this  center 
are  outlined. 
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SUMMARY 


The  research  supported  by  this  grant  is  directed  toward  too. 
general  problems  of  detection  and  identification  of  underground 
explosions  through  the  study  of  radiated  seismic  waves.  Particular 
emphasis  is  on  the  collection  and  analysis  of  broadband  seismic  data 
aL  near  and  regional  distances.  Specific  elements  of  the  research 
program  are:  1)  recording  of  broadband  data  from  events  at  the  Nevada 
Test  Site;  2)  analysis  of  the  coherence  of  ground  motion  near  explosions 
and  earthquakes;  3)  study  of  the  relative  isotropic  and  non-isotropic 
components  of  explosive  sources  through  the  application  of  moment  tensor 
inversion  techniques;  4)  analysis  of  regional  surface  wave  data  in  order 
to  obtain  models  for  the  velocity  and  attenuation  in  the  crust;  5) 
archival  of  near  and  regional  data  sets  which  are  of  value  to  the  general 
discrimination  problem. 

Research  on  elements  1,  3,  and  4  above  is  described  in  the  technical 
report  for  the  first  year  of  this  grant.  Some  of  these  results  have  already 
been  submitted  for  publication  and  should  appear  soon. 

Research  on  element  2,  the  analysis  of  coherence  of  ground  motion, 
is  described  in  sections  II,  III,  IV,  and  V  of  this  report.  This  work 
is  primarily  the  work  of  Keith  McLaughlin  and  forms  part  of  his  PhD 


dissertation. 
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Section  II  describes  the  analysis  of  array  data  recorded  1.9  km 
from  the  explosion  Liptauer  in  Yucca  Valley  of  the  Nevada  Test  Site. 
Coherency  declines  gradually  with  both  inter-station  spacing  and  frequency. 
Broadband  correlation  across  the  400  m  array  is  greater  than  78%  for  the 
horizontal  components  and  greater  that  64%  for  the  vertical  components. 

Section  III  presents  a  new  way  of  looking  at  frequency-wavenumber 
spectral  estimation  with  array  data.  The  differences  between  the 
conventional  beamforming  method  of  estimation  and  the  high-resolution 
method  of  estimation  is  quite  clear  when  expressed  in  terms  of  the 
eigenvalues  of  the  cross  spectrum.  The  conventional  method  uses  primarily 
the  maximum  eigenvalues  while  the  high-resolution  method  uses  primarily 
the  minimum  eigenvalues. 

Section  IV  examines  the  problem  that  exists  when  the  velocity 
structure  is  complicated  to  the  extent  that  it  can  be  considered  to  have 
a  random  component.  This  gives  rise  to  a  random  component  on  the 
seismogram  which  is  usually  referred  to  as  signal-generated  noise  or 
coda.  This  random  component  can  be  separated  from  the  deterministic 
part  of  the  seismogram  and  treated  as  the  convolution  of  the  source  with 
a  stochastic  Green's  function.  Estimates  of  this  stochastic  Green's 
function  and  its  variance  are  useful  in  inferring  bias  and  uncertainty 
in  estimated  source  functions. 

Section  V  treats  the  problem  of  scattering  of  elastic  waves  by 
small  inhomogeneities.  The  solution  is  expressed  in  terms  of  a  moment 
tensor  expansion  of  the  properties  of  the  scatterer.  This  approach  is 
convenient  for  examining  the  trade-offs  between  shape,  heterogeneity, 
and  anisotropy  of  scatterers.  It  reveals  that  a  general  scatterer 
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can  not  be  modeled  by  a  homogeneous  scatterer. 

Progress  on  element  5  of  the  research  program  is  described  in  section 
VI.  In  order  to  provide  better  access  to  archived  seismic  data  and 
facilitate  computations  with  these  data,  a  Computational  Center  for 
Seismology  (CCS)  has  been  established.  The  organization  and  initial 
efforts  of  this  center  are  outlined. 


ANALYSIS  OF  ARRAY  DATA  FOR  EXPLOSION 


LIPTAUER  IN  YUCCA  VALLEY 


ABSTRACT 

The  explosion  Liptauer  (ML  4.7,  BRK)  in  Yucca  Valley,  was  investigated  with  a  small 
array  of  accelerometers  at  an  epicentral  distance  of  5  source  depths  (1.89  km).  The  Yucca  Val¬ 
ley  site  has  pronounced  high  velocity  basement  relief.  A  possible  significant  basement  ofTset 
lay  between  the  array  and  source.  Wavenumber  spectra,  broadband  cross-correlation,  bandpass 
cross-correlation,  and  particle  motion  plots  were  used  to  explore  the  nature  of  the  wave  propa¬ 
gation. 

The  apparent  velocity  of  the  initial  P*wave  at  the  array  was  very  high  (exceeding  20 
km/sec)  for  a  distance  of  1.9  km.  Later  arrivals  on  the  vertical  component  show  a  lower  velo¬ 
city  of  1.2  km/sec.  The  S  waves  at  this  site  and  distance  exhibit  complicated  behavior  Three 
separate  apparent  S  wave  arrivals  with  horizontal  particle  motions  of  distinct  SV,  distinct  SH, 
and  mixed  SH-SV  rectilinear  are  observed.  The  slowest  arrivals  however,  show  no  evidence  of 
lateral  refraction.  Resolution  of  the  arrival  azimuth  for  the  faster  waves  is  insufficient  to  rule 
out  lateral  refractions  as  an  explanation  of  the  transverse  motions.  A  deviatoric  source  as  well 
as  conversions  near  the  source  or  at  dipping  interfaces  are  likely  causes  of  the  strong  transverse 
horizontal  signals. 

The  three  components  of  motion  can  be  ranked  as  radial,  transverse,  and  vertical  increas¬ 
ingly  incoherent.  The  decay  of  inter-station  coherency  with  increasing  station  aeperaiion  is 
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most  pronounced  upon  the  vertical  component.  Broadband  correlation  of  the  radial  and 
transverse  components  was  at  a  minimum  of  78%  across  the  400  meter  array.  The  vertical 
component  of  motion  could  be  found  to  reach  a  minimum  correlation  of  64%.  Bandpass 
filtered  cross-correlations  of  the  vertical  acceleration  components  show  a  steady  decline  in  the 
interstation  coherency  with  increasing  frequency.  No  pronounced  frequency  cutoff  is  evident 
with  the  frequency 'dependence  of  the  interstation  cross-correlation.  The  vertical  acceleration 
signal-to-noise  ratios  for  the  2.5  second  seismogram  are  less  than  than  2-to-l  at  frequencies 
above  10  Hz 

Use  of  the  time  variance  of  the  squared  modulus  (VSM)  was  explored  as  a  measure  of 
the  spread  of  the  seismic  cross-correlation  functions.  This  measure  of  the  cross-correlation 
functions  was  found  to  be  nearly  independent  of  the  bandwidth  up  to  20  Hz  and  only  weakly 
dependent  on  the  cross-correlation  maxima. 

THE  LOCALITY 

Liptauer  was  an  ML  4.7  (BRK),  (ISC  4.8)  explosion  in  Yucca  Valley 

(37.147M16.082*.  April  3,1980,14:00:00.1  UT,  surface  elevation  1335  m.  depth  417  m). 
The  array  of  9  stations  was  located  roughly  5  source  depths  from  the  event  (1.89  km  from  sta¬ 
tion  1,  see  Figure  1).  The  local  geology  is  depicted  in  Figure  2  based  on  reports  by  Barnes 
et  al.  (1963),  and  Colton  and  McKay  (1966).  The  local  stratigraphy  consists  of  a  layer  of  allu¬ 
vium  over  Tertiary  tuffs  filling  a  fault-controlled  valley  of  Paleozoic  sedimentary  basement 
rocks.  The  depth  to  Paleozoic  basement  is  variable  and  probably  only  150  meters  beneath  the 
array.  The  tuffs  may  not  be  represented  directly  beneath  the  array.  The  Carpetbag  Fault  scarp 
is  projected  midway  between  the  array  and  the  source.  The  probable  location  of  the  carpetbag 
fault  can  be  located  by  a  gravity  gradient  in  the  area  where  the  local  Paleozoic  basement 
deepens  sharply  to  the  east.  P-wave  velocities  of  typical  tuff  in  the  area  are  strongly  dependent 
upon  water  content  and  porosity  but  expected  to  be  between  2  and  3  km/ sec  near  shot  depth 
(Keller,  I960).  The  Paleozoic  limestones  are  expected  to  have  P-wave  velocities  exceeding  4.5 
km/ sac. 
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THE  EXPERIMENT 

The  array  of  nine  three-component  force-balance-servo-accelerometers  was  arranged  in  a 
two-dimensional  pattern  of  nested  triangles  (Figure  3).  One  station  failed  to  record  and  is  not 
shown  in  Figure  3.  Horizontal  components  were  aligned  radial  and  transverse  to  the  shot 
azimuth  as  in  Figure  3.  Each  sution  recorded  accelerations  at  200  samples/sec/channel  Five 
pole  anti-alias  filters  were  operated  with  corner  frequencies  at  either  25  or  5f  •_*  indivi¬ 
dual  event  recorders  were  triggered  by  a  common  signal  and  common  relative  avail¬ 

able.  In  Figure  3,  the  relative  arrival  times  (with  respect  to  sution  1)  of  the  P  wave  are 
denoted  in  parentheses  for  each  sution.  Average  P-wave  travel  time  to  the  arty  vas  0  73 
seconds.  The  slowness  of  the  initial  P  wave  break  was  less  than  0.05  sec/km.  Since  such  a 
small  slowness  exceeds  the  theoretical  resolution  of  the  array  for  frequencies  less  than  20  Hz, 
the  traces  were  aligned  on  the  P  wave  break  to  remove  relative  local  sution  delays  for  subse¬ 
quent  analysis 

The  peak  ground  accelerations  (PGA)  for  each  sution  are  ubulated  in  Table  1  and  the 
vertical  PGA’s  are  denoted  by  each  sution  in  brackets  in  Figure  3.  The  scatter  of  the  PGA's 
reflect  the  variation  that  may  occur  in  acceleration  records  over  very  short  disunces.  In  gen¬ 
eral,  the  transverse  PGA’s  are  indistinguishable  from  the  radial  PGA  values.  The  transverse 
peaks  all  come  from  the  same  coherent  transverse  pulse  The  radial  PGA’s  occur  at  different 
peaks  in  the  records. 

Figure  4a,b,c  shows  the  acceleration  records  for  the  array.  Vertical  records  are  relatively 
simple  compared  to  the  horizonul  records.  The  transverse  and  radial  components  show  motion 
coincident  with  the  vertical  first  motion  and  each  sution  is  consistently  positive  radial  and  nega¬ 
tive  transverse.  The  transverse  first  motion  appears  to  have  about  1/2  of  the  predorr  ns.n 
period  of  the  radial  first  motion  and  about  the  same  amplitude. 

The  largest  vertical  amplitude  arise*  from  an  arrival  0.4  second  following  the  P  wm  won 
a  period  of  1/2  second  This  an  va!  is  wet!  developed  on  the  radial  (radial  away*  common* nt 
and  less  prominent  on  the  transverse  (negative  transverse).  Z  versus  R  and  R  veii«s  T 
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acceleration  particle  motions  are  plotted  in  Figures  5A,B,C,  and  D.  Examination  of  the  Z,R 
panicle  motion  reveals  that  the  P-wave  acceleration  is  retrograde  while  the  horizontal  (R,T)  P- 
wave  particle  motion  has  a  slight  counterclockwise  sense. 

A  significant  signal  arrives  upon  the  radial  component  1.0  second  after  the  P  wave  fol¬ 
lowed  by  a' large- transverse  pulse  arriving  1.25  second  after  the  P  wave.  These  motions  are 
clearly  seen  in  Figure  5.  Beginning  1.5  seconds  after  the  P  wave,  the  radial  and  transverse 
accelerations  are  well  correlated,  with  correlation  of  the  +R  and  -T  directions.  The  third 
arrival  possesses  rectilinear  polarization  oriented  nearly  45°  to  the  R  and  T  components. 
These  three  separate  intervals  of  horizontal  particle  accelerations  are  evident  in  the  0.5  secc 
intervals  labeled  1.5  and  2.0  seconds  in  Figure  5. 

If  these  three  separate  arrivals  are  interpreted  as  S  waves  then  it  is  necessary  to  exp) 
the  0.25  second  delay  between  the  apparent  SV  and  SH  waves  as  well  as  the  delay  between  t>.- 
SH  and  rectilinear  SH  and  SV  motion  of  0.25  second.  If  SH-SV  velocity  anisotropy  is  responsi¬ 
ble  for  either  delay,  then  the  velocities  must  differ  by  15%  over  2  km.  P-to-SV  conversion  near 
the  source  could  be  responsible  for  the  SV-SH  delay  as  well  as  a  source  function  with  delay 
between  deviatoric  and  explosive  parts.  The  horizontal  P-wave  motions  are  at  most  1 5  degrees 
off-azimuth  as  inferred  from  the  first  0.5  second  of  R  and  T  motion  in  Figures  SC,  and  D.  If  P 
and  S  ray  paths  are  similar,  it  seems  unlikely  that  the  transverse  component  is  an  off-azimuth 
S-wave  arrival  with  nearly  perfect  transverse  motion.  Furthermore,  the  source  and  receiver  lie 
across  the  strike  of  the  predominant  structure  and  lateral  refractions  should  be  nil  if  the 
geometry  is  truly  2-dimensional. 

Representative  acceleration  spectral  amplitudes  for  a  5.12  second  window  of  all  three  com¬ 
ponents  are  shown  in  Figures  6a  through  6i.  Vertical  and  radial  acceleration  spectra  are  peaked 
between  1.5  and  2.0  Hz,  with  roughly  w*  slopes  below  1.5  Hz.  The  transverse  acceleration 
spectra  have  roughly  a  slope  below  l.S  Hz.  Consequently,  the  transverse  displacement  spec¬ 
tra  have  a  1  slope  and  do  not  have  a  well  determined  low  frequency  asymptote.  This  is  in  con¬ 
trast  to  the  vertical  and  radial  displacement  spectra  which  are  nearly  fiat  below  1.5  Hz.  A  com- 
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1* 

posite  plot  of  transverse  acceleration  spectral  amplitudes  is  shown  in  Figure  7  Slopes  or  2  0 
and  3  0  are  shown  on  the  plot  for  comparison. 

The  ratios  of  the  acceleration  amplitude  spectral  peaks  for  this  window  are  R/T-1.3  +/- 
0  2  and  T/Z-2.2+/-0.1.  The  transverse  signal  is  nearly  as  large  as  the  radial  signal  while  the 
total  vertical  acceleration  signal  is  roughly  1/2  of  the  horizontal  signals.  Spectral  amplitude 
variations  between  stations  are  not  pronounced  below  10  Hz  and  the  spectra  of  one  station 
easily  overlay  the  spectra  of  another.  For  illustration,  the  vertical  spectral  amo.i  -d;*  'ano's  of 
three  stations  are  shown  in  Figure  8a,b,c.  The  amplitude  ratios  are  nearly  flat  from  1  to  10  Hz 
and  become  gradually  more  erratic  with  increasing  frequency. 

F-K  ANALYSIS 

High  resolution  (HR)  frequency  wavenumber  (f-k)  power  spectral  estimates  (Capon 
tt  ah  1969)  were  made  for  a  2  56  second  window  on  all  three  components  of  motion  at 
selected  frequencies.  The  2.56  second  window  encloses  nearly  all  the  significant  signal  beginmg 
with  the  P  arrival.  The  impulse  response  of  the  array  is  seen  in  Figure  9  While  the  main  lobe 
of  the  array  measures  about  1  cycle/km  wide,  the  aliasing  wavenumber  is  ibout  6  cycles/km 
The  sidelobe  pattern  for  this  sparse  array  is  particularly  troublesome.  The  missing  element  of 
the  array  degrades  resolution  in  the  source  direction  and  produces  fou;  protrusions  on  the  main 
lobe  HR  f-k  estimates  have  the  advantage  over  conventional  f-k  estimates  of  suppressing 
some  of  these  features  of  the  impulse  response. 

HR  f-k  power  spectra  for  2.56  seconds  of  the  vertical,  radial,  and  transverse  acceleration 
are  shown  in  Figures  10,11,  and  12  at  3.2,  4.0  and  5  6  Hz  respectively  All  plots  are  the  same 
wavenumber  scale  and  vary  in  slowness  resolution  directly  proportional  to  the  frequency  The  1 
sec/km  slowness  circle  is  labeled  at  each  frequency.  Because  the  records  have  all  been  shifted 
to  align  the  P-wave  arrival,  the  center  of  each  plot  corresponds  to  the  slowness  of  the  P  wave. 
The  P  wave  had  a  slowness  less  than  0.05  sec/km,  therefore  the  mi^location  of  the  origin  of 
each  plot  is  at  most  0.05  sec/km  The  convention  used  is  for  signal  energy  contours  to  plot  at 
the  azimuth  from  which  the  waves  came.  The  azimuth  of  the  shot  is 


indicated  by  the  arrow  on  each  plot. 

The  vertical  component  f-k  spectral  estimate  for  3.2  Hz  show  an  elongation  toward  the 
source  azimuth,  out  to  l.S  sec/km.  Such  large  slownesses  would  be  aliased  at  5.6  Hz  and 
indeed  begin  to  show  a  wrap  around  effect  at  5.6  Hz.  The  vertical  acceleration  records  evi¬ 
dently  contain  a  very  slow  contingent  of  on-azimuth  arrivals  as  well  as  faster  waves  plotting 
near  the  origin.  The  loss  of  contrast  of  the  vertical  spectral  peak  above  the  background  from 
14db  at  3.2  Hz  to  8  db  at  S.6  Hz  continues  at  higher  frequencies.  The  radial  f-k  spectra  is  com¬ 
pact  at  3.2  Hz  and  mimics  the  impulse  response  with  a  contrast  of  28  db  above  the  background. 
At  4.0  and  S.6  Hz,  the  radial  spectra  are  elongated  along  the  source  azimuth  but  are  dominated 
by  the  high  apparent  velocity  energy  that  plots  near  the  origin.  The  radial  component  may  con¬ 
tain  waves  as  slow  as  0.5  sec/km  with  spectral  power  12  db  down  from  the  much  faster  compli¬ 
ment  of  signal.  The  transverse  spectra  are  elongated  in  the  direction  of  the  source  at  3.2  Hz 
but  compact  and  near  the  origin  at  4.0  to  5.6  Hz  and  at  other  intermediate  frequencies  not 
shown.  The  spreads  of  the  f-k  spectra  perpendicular  to  the  source  azimuth  are  clearly  limited 
by  the  width  of  the  central  lobe  of  the  impulse  response  shown  in  Figure  9.  If  we  interpret  the 
f-k  power  spectra  as  the  Fourier  transform  of  the  spatial  covariance  function  then  this  means 
that  the  spatial  correlation  function  is  much  broader  than  the  array  dimension  at  these  frequen¬ 
cies  The  vertical  component  f-k  spectral  maximum  falls  from  14  db  above  the  background  at 
3.2  Hz  to  only  9  db  above  the  background  at  S.6  Hz.  This  trend  continues  with  increasing  fre¬ 
quency  implying  that  the  vertical  acceleration  field  has  declining  coherency  with  increasing  fre¬ 
quency. 

INTER-STATION  BROADBAND  CORRELATIONS 

It  is  necessary  to  explore  the  coherency  of  the  seismic  traces  without  applying  phase  shifts 
corresponding  to  specific  slownesses  and  azimuths  required  to  align  the  records,  because  the 
record  windows  used  enclose  more  than  a  single  arrival  with  well  determined  slowness.  To 
accomplish  this  broadband  acceleration  cross-correlations  were  computed  for  the  same  2.S 
second  window  used  in  the  f-k  analysis.  Examples  of  the  cross  correlation  functions  can  be 
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seen  in  Figures  13a, b,c  The  maxima  of  the  cross-correlations  are  plotted  in  Figures  14a, b,c 
and  ISa.b.c  The  variance  of  the  squared  modulus  (VSM)  is  also  listed  in  table  2  as  a  measure 
of  the  width  of  each  correlation  function  A  discussion  of  the  VSM  is  given  in  Appendix  C. 
Best  lapse  rates  (distance  at  which  the  correlation  falls  to  1/e)  for  the  transverse  and  radial 
directions  were  fit  to  the  data  and  the  correlation  contours  are  plotted  for  comparison  with  the 
data  in  Figures  ISa.b.c.  The  lapse  rates  across  the  wave  front  were  2600  m,  6800  m,  and  S200 
m  for  the  vertical,  radial,  and  transverse  components  respectively.  In  the  radial  direction,  the 
lapse  rates  were  672  m,  4000  m,  and  1600  m  for  the  vertical,  radial  and  transverse  components. 
In  contrast  to  the  Colwick  experiment,  the  Liptauer  array  data  are  very  coherent  over  the  400 
meter  array.  The  minimum  correlations  for  the  three  components  are  0.64,  0.84,  and  0.78  for 
the  vertical,  radial,  and  transverse  components,  respectively.  The  decay  of  correlation  with  sta¬ 
tion  separation  is  much  less  pronounced  than  was  observed  for  the  Colwick  array  where  a  simi¬ 
lar  time  window  might  exhibit  correlations  of  0.S  for  stations  separated  only  400  meters.  The 
two-dimensional  character  of  the  correlation  is  strongly  developed,  showing  much  different 
decay  rates  for  the  two  orthogonal  directions. 

The  inter-station  correlations  for  different  components  can  be  seen  in  the  scattergrams  of 
R*  maxima  versus  Rr  maxima  and  RT  maxima  versus  the  Rz  maxima  (Figure  16).  For  any 
given  inter-station  pair  ,  radial  components  are  better  correlated  than  the  transverse  com¬ 
ponents,  and  transverse  components  are  better  correlated  than  the  vertical  components  R  is 
more  coherent  than  T,  and  T  is  more  coherent  than  Z,  while  the  total  R  signal  strength  is 
greater  than  T,  and  similarly  the  total  T  signal  strength  is  larger  than  Z.  Ratios  of  the  total  sig¬ 
nal  strengths,  as  measured  by  the  spectral  peaks  at  1.5  Hz,  are  R/T  —1.3  and  T/Z  —2.2.  This 
consistency  between  components  suggests  that  lack  of  coherency  is  due  to  a  common  noise  sig¬ 
nal  on  all  three  components.  The  noise  signal  would  become  less  correlated  with  increasing 
separation  between  station  pairs  and  give  predictable  ratios  of  relative  aignal-to- noise  between 
the  three  components.  We  can  define  signal- to- noise  power  ratios  from  the  correlation  func¬ 


tions, 
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^I+lNj/Sr)1  **“l+(Az/Sz)2 

where  Sn/N *,  Sr/ST,  end  Sz/Nz  we  the  radial,  transverie,  and  vertical  average  signal-to-noise 
ratios.  The  'noise*  is  by  definition  the  signal  component  not  common  to  the  two  signals  being 
correlated  In  Figure  17  the  inferred  values  of  these  signal-to-noise  ratios  are  plotted  for  radial 
versus  transverse  and  transverse  versus  vertical  for  each  station  pair.  If  we  assume  the  the 
errors  in  the  inferred  signal-to-noise  ratios  we  distributed  evenly  for  radial,  transverse,  and 
vertical  signals,  then  ratios  of  ( N*IS* ) V ( Nt/St) j.  and  (A,r/5j-),/(Arz/£z)2  we  (0.8 1 ) 3  and 
(0.64) 2  respectively.  This  is  consistent  with  the  spectral  ratios  of  R/T  and  T/Z  of  1.3  and  2.2 
near  1.5  Hz.  The  R  and  T  components  have  newly  equal  noise  component,  and  the  vertical 
component  has  roughly  1/2  as  much  noise  power  as  the  horizontal  components. 

FREQUENCY  DEPENDENCE  OF  SPATIAL  CORRELATION 

Bandpass  filtering  of  the  vertical  acceleration  cross-correlation  functions  was  performed 
for  frequency  bands  of  1.25-2.5,  2.5-5.0,  5.0-10.0,  and  10.0-20.0  Hz.  The  maxima  of  these 
cross-correlation  functions  are  listed  in  Table  3.  Examples  of  the  bandpass  filtered  auto-  and 
cross-coirelation  functions  we  shown  in  Figure  19A,B  for  stations  8  and  9  and  1  and  7.  The 
cross-correlation  functions  of  stations  8  and  9  we  more  symmetric  than  the  correlation  func¬ 
tions  of  stations  1  and  7.  Speciflcly,  the  5.0-10.0  Hz  cross-correlation  function  for  stations  1 
and  7  (Figure  19B)  is  biased  toward  positive  time  delays.  Stations  1  and  7  have  the  largest 
range  sepwation  of  any  station  pair  for  the  array.  This  asymmetry  presumably  corresponds  to 
the  propagation  of  the  vertical  waveforms.  This  asymmetry  is  not  as  prominent  at  the  lower 
frequencies.  Attempts  to  quantify  such  an  observation  with  statistics  such  as  the  centroid  of 
the  modulus  or  square  modulus  of  the  correlation  function,  or  the  location  of  the  maximum 
peak  were  not  fruitful.  One  reason  for  this  is  the  uncertainty  of  the  correlation  peak,  or  peaks, 
expreased  by  the  variance  of  the  squared  modulus  is  dismissed  in  the  Appendix. 

The  maxima  of  these  bandpass  filtered  correlation  functions  are  plotted  in  Figure  19. 
Maximum  correlations  are  all  above  0.75  for  bandwidths  below  2.5  Hz  and  above  0.5  for 
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bandwidths  below  5.0  Hz.  The  data  favor  a  gradual  decline  of  correlation  with  frequency  rather 
than  an  abrupt  decline  as  seen  at  the  Colwick  array.  All  correlations  have  fallen  to  0.67  or  less 
above  10  Hz.  This  would  correspond  to  a  signal-to- noise  ratio  below  2-to-l  at  this  bandwidth  of 
10  to  20  Hz 

DISCUSSION 

The  lack  of  simple  'layer-cake'  structure  in  Yucca  Valley  surely  has  profound  implications 
for  use  of  acceleration  records  such  as  those  at  the  Liptauer  array.  Of  ^articular  interest  is  the 
origin  of  the  horizontal  accelerations  observed  at  1.89  km  from  an  mt  4.7  explosion.  It  is  not 
clear  from  the  data  available,  that  the  coherent  transverse  signals  require  a  deviatoric  source. 
The  transverse  component  is  very  coherent  across  the  array,  less  than  1/4  of  the  2.S  second 
seismogram  can  be  interpreted  as  noise.  The  predominant  transverse  signal  begins  with 
apparent  SH  motion  followed  by  particle  motion  in  phase  with  the  radial  component.  The 
transverse  slowness  spectra  below  6  Hz  is  consistent  with  on-azimuth  arrivals  at  high  apparent 
velocities.  Only  the  initial  transverse  motion  can  be  explained  as  P-SV  motion  at  10  to  20 
degrees  off-azimuth  with  the  source. 

The  simplest  explanation  for  the  transverse  pulse  is  SH  wave  generation  near  the  source 
arriving  at  the  array  with  a  high  apparent  velocity.  The  low  frequency  amplitude  spectra  of  the 
transverse  signal  differ  from  the  radial  and  vertical  spectra  (Figures  6  and  7).  The  transverse 
component  has  a  disproportionate  lack  of  low  frequencies  below  1.5  Hz  compared  to  the  radial 
and  vertical  components.  With  uncertainties  in  the  local  propagation  it  is  not  possible  to 
address  the  possibility  that  the  SH  wave  generation  has  a  different  frequency  dependence  than 
the  P-SV  wave  generation.  If  the  SH-SV  delay  of  0.25  seconds  were  taken  for  face  value,  it 
might  be  argued  that  the  source  of  the  SH  waves  is  delayed  with  respect  to  the  P-SV  source. 
However,  it  seems  likely  that  the  initial  radial  signal  is  a  P-to-SV  conversion  near  the  source. 
The  remaining  rectilinear  S-wave  motion  polarized  45°  to  the  R  and  T  directions  and  delayed 
0.25  sec  to  the  SH  arrval  complicates  the  picture.  The  avalible  data  can  not  establish  whether 
this  arrival  represents  S-wave  motion  with  a  separate  ray  path,  anisotropic  motion,  or 


13 


conversion  in  the  high  velocity  basement.  F-k  analysis  shows  that  slow  arrivals  are  on  azimuth, 
but  does  not  have  sufficient  resolution  to  rule  out  steeply  incident  off-azimuth  S  waves.  The 
high  apparent  velocities  of  the  radial  and  transverse  components  (Figures  10,11,  and  12)  as  well 
as  the  small  vertical  amplitudes  (Figure  4)  indicate  that  steeply  incident  S  waves  from  the  base¬ 
ment  are  important.  The  structural  inhomogeneity  of  the  Paleozoic/Tertiary  contact  beneath 
Yucca  Valley  (Figure  2)  makes  any  of  these  interpretations  uncertain.  The  deployment  of  an 
accelerometer  array  in  Yucca  Valley  has  shown  that  the  near-source  seismogram  can  possess  a 
wealth  of  complexity. 

Broadband  signals  at  the  Liptauer  array  were  very  coherent  and  consistent  with  the  model 
for  a  common  noise  signal  of  equal  size  on  the  radial  and  transverse  components  while  1/2  as 
large  on  the  vertical  component.  This  noise  signal  must  have  a  correlation  length  exceeding  the 
array  dimension  of  400  meters  up  to  6  Hz,  otherwise  the  spatial  cross-correlations  would  reach 
an  asymptotic  level  within  the  array.  Spatial  correlation  declines  steadily  with  increasing  fre¬ 
quency.  The  spatial  and  frequency  dependence  of  the  vertical  cross-correlations  have  been 
combined  in  Figure  20A.  Some  suggested  contours  are  plotted.  For  contrast,  a  similar  plot  of 
the  Colwick  data  is  shown  in  Figure  20B.  The  primary  difference  between  the  two  cases  can  be 
seen  at  the  lower  frequencies.  The  decline  of  coherency  at  the  Colwick  array  had  a  sudden 
onset  near  5  Hz,  while  the  Liptauer  data  favors  a  more  gradual  decline  in  coherency  with 
increasing  frequency. 
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TABLE  1:  UPTAUER  PGA’l  (g) 

STA  Z  B  T 

1  0.092  0.067  0.093 

2  0.060  0.067  0.100 

3  0.010  0.093  0.017 

4  0.06S  0.073  0.067 

6  -0.062  -0.060  0.073 

7  0.049  -0.100  0.073 

g  0.065  0.080  0.067 

9  0.060  0.080  0.080 
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TABLE  2:  BROADBAND  CROSS-CORRELATION  MAXIMA  AND  VSM 


pftir 

Z 

R 

T 

Z  (VSM) 

R  (VSM) 

T  (VSM) 

»ep*rttion 

1.2 

.88 

.97 

.93 

.23 

'.24 

.20 

100 

1.3 

.91 

.97 

.97 

.36 

.26 

.21 

too 

1.5 

.81 

.93 

.91 

.32 

.26 

.19 

200 

1.6 

.69 

.95 

.87 

.46 

.26 

.22 

200 

1.7 

.70 

.88 

.78 

.34 

.24 

.20 

346 

1.6 

.83 

.96 

.93 

.27 

.26 

.20 

200 

1.9 

.84 

.95 

.87 

.27 

.26 

.17 

173 

2.3 

.90 

.92 

.96 

.26 

.27 

.20 

100 

2.5 

.92 

.95 

.92 

.18 

.27 

.18 

173 

2.6 

.75 

.98 

.92 

.37 

.28 

.20 

100 

2,7 

.74 

.90 

.84 

.26 

.26 

.18 

265 

2,8 

.86 

.97 

.96 

.16 

.27 

.19 

173 

2.9 

.88 

.98 

.93 

.16 

.27 

.15 

100 

3,5 

.81 

.92 

.89 

.34 

.28 

.19 

265 

3,6 

.64 

.95 

.86 

.55 

.29 

.22 

173 

3,7 

.65 

.84 

.80 

.41 

.29 

.20 

265 

3,8 

.90 

.98 

.97 

.28 

.27 

.20 

100 

3,9 

.85 

.97 

.88 

.30 

.28 

.28 

100 

5,6 

.78 

.96 

.91 

.38 

.29 

.29 

200 

5,7 

.73 

.90 

.82 

.30 

.28 

.16 

400 

5.8 

.74 

.92 

.91 

.21 

.28 

.18 

346 

5,9 

.82 

94 

.89 

.20 

.29 

.13 

265 

6,7 

.95 

.92 

.93 

.27 

.28 

.18 

200 

6,8 

.67 

.95 

.93 

.32 

.29 

.19 

200 

6.9 

.85 

.98 

.97 

.29 

.30 

.16 

100 

7.8 

.65 

.84 

.88 

.24 

.28 

.17 

200 

7,9 

.81 

.88 

.93 

.22 

.29 

.14 

173 

8,9 

.91 

.98 

.94 

.17 

.28 

.15 

100 
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TABLE  3:  BANDPASS  CROSS-CORRELATION  MAXIMA  AND  VSM  (jrc1) 


station 

1.25- 

2.5- 

5.0- 

10.0- 

broad¬ 

d 

pair 

1.3 

2.5 

5.0 

10.0 

20.0 

band 

(m) 

.93 

.90 

.76 

.33 

.91 

100 

.39 

.34 

.. 

.56 

.36 

200 

1.5 

.15 

.71 

.59 

.43 

.82 

.34 

.32 

.56 

— 

.32 

200 

1,6 

.97 

.55 

.55 

.42 

.69 

.45 

.55 

.32 

.43 

.46 

1.7 

.15 

.54 

.52 

.52 

.43 

346 

.35 

.40 

.29 

.34 

.34 

264 

3,7 

.76 

.53 

.53 

.34 

.65 

.42 

.43 

— 

.44 

.41 

200 

5,6 

1.00 

.68 

.74 

.66 

.78 

.37 

.53 

-- 

.46 

.38 

100 

*.9 

.96 

.80 

.79 

.67 

.91 

.16 

.27 

.38 

.51 

.17 
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FIGURE  CAPTIONS 

Figure  1.  Location  of  Liptauer  and  the  accelerometer  array.  The  Carpetbag  Fault  trace  is  indi¬ 
cated. 

Figure  2.  Diagrammatic  cross-section  of  the  local  geology  through  the  array  and  shot  point  of 
Figure  1. 

Figure  3.  Array  geometry,  definition  of  radial  (R)  and  transverse  (T)  motions,  and  direction  of 
the  shot  are  indicated.  Relative  P  times  and  vertical  peak  accelerations  in  g*s  are 
annotated  in  parentheses  and  brackets  respectively. 

Figure  4A.B.C.  Vertical  (A),  radial  (B),  and  transverse  (C)  acceleration  traces  recorded  at  the 
array. 

Figure  SA,B,C,D.  Acceleration  particle  motions  for  (A,B)  Z  versus  R,  and  (C,D)  R  versus  T 
for  the  8  stations  of  the  array.  The  Z  versus  R  plots  have  an  exaggeration  of  50%  in 
the  R  direction  relative  to  the  Z  direction,  the  V  marks  the  beginning  of  the 
record.  Each  0.5  seconds  of  record  is  portrayed  separately.  Records  do  not  start  at 
the  origin  due  to  small  D.C.  offsets  in  each  recording. 

Figure  6A-I.  Log-log  acceleration  spectral  amplitude  plots  for  the  2.56  seconds  of  data  on  each 
component  of  motion. 

Figure  7.  A  composite  plot  for  the  8  transverse  acceleration  spectra.  A  2.0  and  3.0  slope  are 
shown  for  comparison. 

Figure  8.  Spectral  ratios  of  vertical  acceleration  spectra  for  stations  3,5,  and  6. 

Figure  9.  Wavenumber  impulse  response  for  the  8  station  array.  The  main  lobe  is  about  1 
cycle/km  wide  and  the  aliasing  wavenumber  in  most  directions  is  about  6  cycle/km. 

Figure  10.  High  resolution  (HR)  frequency-wavenumber  (f-k)  power  spectral  estimate  for  the 
Z,R,T  components  at  3.2  Hz.  Contours  are  1,2,  or  3  db  with  respect  to  the  max¬ 
imum.  The  1  sec/km  slowness  circle  is  indicated  and  the  shot  azimuth  is  indicated 
by  the  arrow.  The  convention  used  is  for  the  arriving  energy  to  plot  at  the  azimuth 
from  whence  it  caifte.  The  Z,  R,  and  T  maxima  are  14,  28,  and  12  db  above  the 
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background  respectively. 

Figure  11.  HR  f-k  power  spectre  for  the  Z,R,T  components  st  4.0  Hr.  The  Z,  R,  end  T  max¬ 
ima  ere  14, 16,  end  IS  db  above  the  background  respectively. 

Figure  12.  HR  f-k  power  spectre  for  the  Z,R,T  components  at  5.6  Hz.  The  Z,  R,  end  T  max¬ 
ima  ere  9,  16,  end  IS  db  above  the  maxima  respectively. 

Figure  13A.B.C.  Broadband  acceleration  cross-correlation  functions  for  stations  3,5,  and  8  for 
(A)  the  Z  component,  (B)  the  R  component,  and  (C)  the  (T)  components.  The 
maxima  of  the  normalized  correlation  functions  are  indicated. 

Figure  14A.B.C.  Cross-correlation  maxima  plotted  versus  inter-station  separation.  Best  fit 
radial  and  transverse  exponential  decay  curves  are  plotted  for  comparison. 

Figure  15A,B,C.  Cross-correlation  maxima  plotted  versus  transverse  and  radial  inter-station 
separation.  Best  fit  exponential  surfaces  are  contoured  for  comparison. 

Figure  16A.B.  A.)  Radial  versus  transverse  correlation  maxima  and  B.)  Transverse  versus 
Vertical  correlation  maxima.  Each  point  represents  a  station  pair. 

Figure  17A.B.  Inferred  signal-to-noise  ratios  for  the  A.)  radial  versus  transverse  and  B.) 

transverse  versus  vertical  components.  Each  point  represents  a  station  pair  from 
Figure  16.  The  best  fit  slope  is  indicated.  The  "noise"  is  defined  as  the  uncommon 
signal  between  the  two  stations. 

Figure  18.  Maxima  of  the  bandpass  filtered  cross-correlation  functions  for  vertical  acceleration 
2. 56  second  window  from  Table  3.  Station  separations  are  indicated  next  to  the  plot¬ 
ted  points  in  meters. 

Figure  I9A.B.  Bandpass  auto-  and  cross-correlation  functions  for  stations  A.)  I  and  7,  and  B.) 

8  and  9.  The  normalized  maxima  of  the  cross-correlation  are  annotated  and  the 
inferred  width  estimated  from  the  variance  of  the  squared  modulus  (VSM)  is  shown 
as  a  horizontal  bar  and  labeled  in  units  of  seconds. 

Figure  20A,B.  A.)  Maxima  of  the  Liptauer  array  vertical  acceleration  cross  correlation  functions 
contoured  on  a  plot  of  inter-station  separation  and  frequency  bandwidth.  B.)  A 
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comparison  plot  of  the  cross-correlation  maxima  for  the  Colwick  array  data. 
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APPENDIX 


VARIANCE  OF  THE  SQUARED  MODULUS  OF  THE 
CROSS-CORRELATION  FUNCTION  AS  A  MEASURE  OF 

THE  CROSS-CORRELATION  FUNCTION 


The  width  of  a  cross-correlation  peak  is  dominated  by  the  cross-spectral  bandwidth  and 
may  not  represent  the  true  uncertainty  of  the  location  of  that  peak.  The  half  width  at  half  max¬ 
imum  (HWHM)  is  often  cited  as  the  uncertainty  of  the  optimum  cross-correlation.  Because  the 
HWHM  represents  the  bandwidth  of  the  cross-spectrum  and  not  the  net  width  of  the  correla¬ 
tion  function,  another  measure  is  desired.  For  example,  two  adjacent  peaks  of  nearly  equal 
height  may  be  separated  by  a  deep  negative  minimum,  and  the  widths  of  the  peaks  do  not 
speak  for  the  uncertainty  implied  by  their  mutual  adjacency.  Bracewell  (1978)  suggests  that  the 
mean  square  departure  from  the  centroid  as  a  measure  of  the  spread  of  a  correlation  function 
with  zero  mean.  The  variance  of  the  squared  modulus  (VSM)  is  defined  as 

where  R(t)  is  the  correlation  function  and  J r|/?  (r)|Jdt  is  the  centroid  of  the  squared  modulus. 
The  VSM  is  tabulated  with  the  correlation  maxima  in  Tables  2,  and  3.  The  correlation  of  the 
VSM  with  the  correlation  maxima  for  the  vertical  component  acceleration  is  shown  in  the  scat- 
tergram  of  Figure  A1  and  with  the  station  separation  in  Figure  A2.  No  separate  measure  of  the 
quality  of  the  correlations  is  available,  but  if  the  VSM  was  to  be  a  useful  statistic  we  would 
expect  it  to  show  the  same  systematic  as  the  correlation  maxima.  The  VSM  are  nearly  con- 


67 


siant  for  a  given  component  and  independent  of  the  bandwidth  as  demonstrated  in  Figures 
'S  ',  and  B  .  These  equivalent  widths  of  the  correlation  functions  are  sh'  r  as 

bars  over  the  maxima  of  the  auto-  and  cross-correlation  functions  in  Figures  19A  and  B 
.  The  use  of  this  measure  of  the  equivalent  width  seems  doubtful  as  a  measure  of 
the  uncertainty  of  the  maxima  locations.  Windowing  of  the  cross-correlation  functions  about 
t>  »ir  centroid  with  a  window  of  the  specified  width  would  however  capture  the  significant  posi¬ 
tive  and  negative  correlations  and  may  be  useful  if  such  an  estimate  of  the  correlation  spread  is 
jc: red 
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FIGURE  CAPTIONS 

Figure  Cl.  Variance  of  the  squared  modulus  (VSM)  versus  maxima  of  the  vertical 
acceleration  cross-correlation  functions  listed  in  Table  2.  Each  point  represents  a  station 
pair. 

Figure  C2.  VSM  versus  distance  of  separation  for  station  pairs  shown  in 
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III 


ANALYSIS  OF  HIGH-RESOLUTION  FREQUENCY- 
WAVENUMBER  SPECTRAL  ESTIMATOR  WITH  NON- 

STATIONARY  SEISMIC  DATA 


Abstract 

Previous  methods  for  analysis  of  frequency-wavenumber  spectral  estimation  have  relied 
on  stationary  and  multidimensional  noise  estimators.  Often,  seismologists  wish  to  estimate 
slowness  spectra  of  pulse-like  waveforms  or  other  types  of  non-sutionary  signals.  The  method 
of  high  resolution  (HR)  frequency- wavenumber  (f-k)  spectral  estimation  is  explored  using 
singular  value  decomposition  of  the  cross-spectrum  estimate.  Analysis  shows  that  the 
difference  between  the  conventional  beamforming  method  and  the  HR  estimator  is  the  weight¬ 
ing  given  the  eigenvalue-eigenvector  contributions  to  the  cross-spectrum.  It  is  shown  that  the 
conventional  method  weights  the  largest  eigenvalues  while  the  HR  method  utilizes  the  smallest 
eigenvalues.  The  near-singular  contributions  to  the  cross-spectral  matrix  give  an  approximate 
view  of  the  null  space  of  the  matrix.  The  HR  estimator  gives  an  inverted  view  of  the  null  space 
of  the  cross-spectral  matrix. 

An  example  from  an  experiment  exhibiting  multipathing  behavior  is  used  to  demonstrate 
the  rank  deficiency  of  real  seismic  phase  delay  data.  The  decomposition  points  out  the  disad¬ 
vantages  and  advantages  to  the  HR  method  for  detection  of  resolvable  multiple  arrivals. 
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Introdectian 

The  high  resolution  (HR)  frequency-wavenumber  (f-k)  spectral  estimator  introduced  by 
Capon,  Greenfield,  and  Kolker  (1967)  was  presented  as  a  minimum  variance,  unbiased, 
maximum-likelihood  filter.  The  filter  passes  undistorted  waveforms  at  the  steering  slowness 
while  optimally  rejecting  noise  power  from  other  slownesses.  Woods  and  Lintz  (1973)  showed 
that  the  increased  resolution  of  the  HR  method  arises  from  the  assumption  of  correlated  plane 
waves.  They  synthesized  2  plane  waves  of  varying  correlation  and  explored  time  window  sam¬ 
pling,  and  the  effects  of  additive  white  noise.  They  demonstrate  that  the  HR  method  can 
resolve  two  closely  spaced  plane  waves  when  conventional  slowness  spectral  estimate  may  not. 
Cox  (1973)  investigated  the  effects  of  noise  and  interfering  arrivals  on  optimal  array  processors, 
and  showed  the  resolving  power  and  effective  gain  of  the  HR  filter  to  be  intimately  related  to 
the  noise  cross-spectrum  and  the  mismatch  of  the  steering  vectors  measured  in  a  metric  defined 
by  the  cross-spectral  matrices  (  noise  and  signal).  The  purpose  of  this  paper  is  to  propose  a 
physical  interpretation  of  the  HR  estimator  without  the  need  for  a  noise  model.  By  the  analyti¬ 
cal  technique  of  singular  value  decomposition  of  the  estimated  cross-spectral  matrix,  the  HR 
method  can  be  seen  to  be  a  best  fit  of  the  phase  delay  data  to  a  superposition  of  plane  waves 
provided  the  interference  of  the  incoming  plane  waves  can  be  reduced  by  averaging.  The  use  of 
eigenvector  (or  principal  component)  decomposition  of  the  cross  spectral  matrix  was  used  by 
Der  and  Flinn  (1975)  where  they  showed  that  two  signals  could  be  independently  resolved  by 
an  array  if  1.)  the  signals  are  not  of  comparable  amplitude  or  2.)  that  each  signal's  slowness 
did  not  coincide  with  the  side  lobes  of  the  beam  directed  at  the  other  signal.  By  examination  of 
the  perturbation  theory  of  matrix  decomposition  we  see  that  some  of  these  peculiarities  of  the 
high  resolution  method  are  related  to  the  perturbations  of  the  cross-spectral  matrix  in  actual 
practice. 

The  time  dependent,  two-dimensional  wave  field,  u(x,y,t),  may  be  represented  by  a  triple 
Fourier  transform  (Burg,  1964) 

«(jr,y,  r)- f  dm J  dkj  dfc,  w ( krm)  e*p\-i(  kxx+ k,y—m  t)  1 


(1) 


t 


? 
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Where  k„  and  k,  are  horizontal  wavenumbers  in  the  x  and  y  directions.  The  corresponding 
inverse  relationship, 

dxj  4yu(xty,u)  exp[i(kx+k,)]  (2) 

is  assumed  to  hold  for  some  region  of  (x,y).  We  observe  the  seismic  field  u(x,y,t)  at  discrete 

points  and  the  frequency-wsvenumber  spectrum  that  duplicates  the  field  is  repetitive 

in  ( kjpk ,)  with  an  aliasing  wavenumber  that  varies  with  azimuth.  The  conventional  method  for 

estimation  of  the  signal  enrgy  spectra,  P( k,<u)  —  |«(k,<u)|J,  is  to  replace  the  spatial  integrals 

with  a  weighted  sum  over  the  sampled  wave  field  after  performing  a  Fourier  transform  upon 

the  sampled  time  domain  window. 

/,c(k,«)-l£ii(ry,M)«apl/kery]|J  (3) 

y-i 

-  £  Uj(u)  1/jW  (fy— f/)l 

/.*- 1 

-  £  S/»)  UjX  k)  (/,<  k) 

-u#(ys(*)u(k) 

where  the  sution  coordinates  of  the  array  are  r„  Suim)  is  the  cross-spectral  matrix  estimate, 
is  the  sution  phase  delay  or  steering  vector,  U*  is  the  Hermitian  transpose 
of  U,  and  u'  is  the  complex  conjugate  of  u.  The  conventional,  or  beamforming  estimate, 
P' r(k,»),  is  the  convolution  of  the  beam  pattern,  /(k),  with  the  true  wavenumber  spectrum, 

«(k.«) 

Pclk,m)»\f  dk7(k')ir(k-k')l1  (4) 

/(k)-^o»t/keryl  (5) 

Kk)  is  the  Fourier  transform  of  the  spatial  sampling  function,  the  response  to  a 

M 

vertically  incident,  k-0,  plana  wave,  u(k)«*(k).  A  beamform,  or  stack  of  the  seismograms, 
from  an  array  is  formed  by  the  sum  over  the  suitably  time  shifted  and  weighted  seismograms. 

**(/■, s)-£»>/f-ser,)  <«> 

The  Fourier  signal  energy  spectrum  of  the  beam  for  slowness,  a,  is  identical  to  the  coovantiooal 
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frequency- wavenumber  spectrum  estimate  for  k-ws,  with  Wj*  1, 

Fc(k,«)«*l3£,  M//u(w)«p[/k*ry]|3  (7) 

7"1 

Consequently,  the  conventional  spectral  estimator  produces  a  superposition  of  shifted  beam  pat¬ 
terns  corresponding  to  plane  waves  required  to  represent  the  field  observed  at  discrete  points, 
the  station  distribution.  The  presence  of  beam  sidelobes  and  multiple  arrivals  can  produce 
complicated  interference  patterns.  A  small  array  with  significant  sidelobes  can  be  difficult  to 
interpret  in  the  presence  of  interfering  arrivals.  To  investigate  the  nature  of  the  estimated  sig¬ 
nal  energy  spectta  requires  calculation  of  Fc(k)  on  a  grid  of  (kx,ky)l  contouring,  and  genera¬ 
tion  of  graphical  output.  Even  in  the  case  of  a  single  arrival  these  calculations  can  be  consider¬ 
able.  A  simple  search  for  the  global  maximum  is  complicated  by  the  many  local  maxima  of  the 
beam  response.  A  primary  advantage  of  the  HR  estimator  is  the  suppression  of  sidelobes  for 
lone  or  non-interfering  multiple  arrivals. 

The  estimation  of  the  cross-spectral  density  matrix,  S0(m)  ■*  u,  («)«/(«),  for  stationary 
noise  data  is  often  performed  by  averaging  over  several  temporal  windows.  This  is  unsatisfac¬ 
tory  for  non-stationary  seismic  data.  An  estimate  for  non-stationary  signals  is  either  made  by 
smoothing  the  cross-spectrum  with  a  convolution  operator,  or  Fourier  transforming  the  win¬ 
dowed  cross-correlation  functions  as  estimates  of  the  spatial  covariance  function.  The  two 
methods  are  equivalent  while  there  may  be  computational  advantages  for  a  narrow  band, 
l/T  -  /.,  frequency  domain  smoothing  operator  over  the  time  domain  windowing  of  the  spatial 
covariance  estimate.  Advantages  include  the  reuse  of  the  cross-spectra  for  coherency  estimates 
and  rewindowing.  In  either  case  the  cross-spectral  matrix  may  be  written  as 

Sy(»  w,u,(m  /)  uj(m  /)  (I) 

A-*, 

where  the  weighted  sum  over  frequencies  wk|  to  mtj  produces  s  smoothed  croes-epectrai  esti¬ 
mate  with  center  frequency 
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The  use  of  such  smoothing  operators  is  straight  forward  and  the  statistics  and  pitfalls  are 
described  in  works  such  as  Jenkins  and  Watts  (1968).  A  narrow  averaging  kernel,  (w/1,  does 
not  stabilize  the  estimate  and  contains  positive  bias  for  the  purpose  of  estimating  coherency 
between  stations.  If  the  spectral  smoothing  is  done  over  too  wide  a  bandwidth,  the  operation 
will  destructively  average  the  slow  deterministic  signals.  An  averaging  kernel  should  be  narrow 
enough  to  admit  any  time  delay  within  interest  across  the  array.  For  a  frequency  smoothing 
operator  of  width  1  Hz,  a  signal  with  slowness  of  2  sec/km  would  be  seriously  degraded  for  sta¬ 
tions  1  km  apart.  The  phase  delay  d“2ir/wi*s  across  the  frequency  band  ,/„,  should  remain 
much  less  than  n.  Further  discussion  of  the  wavenumber  power  spectra  assumes  a  cross- 
spectral  matrix,  Su,  has  been  made. 

The  HR  f-k  estimator 

The  high  resolution  (HR)  estimate  of  the  signal  energy  spectrum  may  be  written 

F**(k>MU"(k')S-,U)U(k')]-‘  (9) 

where  S-1  is  the  inverse  of  the  cross-spectral  matrix.  Capon  et.  al.  (1967)  introduced  the  esti¬ 
mator  as  a  distortionless  filter  for  k-k\  while  optimally  rejecting  signal  power  at  k*k'.  For  a 
multidimensional  Guassian  noise  distribution  the  filter  is  a  maximum  likelihood,  minimum 
variance  estimate.  The  HR  estimate  requires  only  the  additional  calculation  of  the  inverse 
cross-spectral  matrix.  The  similarity  of  equation  9  to  equation  3  is  clear. 

Cox  (1973)  decomposes  Su  into  the  noise  cross-spectral  matrix,  Qv,  and  the  signal  cross- 
spectral  matrix,  Rv,  S-Q4-R.  The  matrices,  S,S~,,R,R'\Q,Q“\  define  metrics  where  any  two 
vectors  a,  and  k  may  be  represented  by  their  expansion  as  eigenvectors.  The  eigenvectors  of 
one  the  matrices  may  not  span  the  entire  space  and  a  portion  of  a,  or  h,  may  not  project  to  the 
eigenvector  expansion.  The  angle  between  a,  and  h  is  defined  by  a  generalized  inner  product 
of  the  eigenvector  expansions  in  these  metrics.  Consequently,  the  ability  to  resolve  any  two 
vectors  is  described  by  the  respective  null  spaces  of  the  cross-spectral  matrices.  Examining  the 
form  of  the  beam  and  HR  estimates  in  equations  3  and  9,  note  that  components  of  the  station 
phase  delay  vectors,  U(k),  belonging  to  the  null  spaces  of  S  and  S'*  make  no  contributions  to 
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the  the  respective  beam  and  HR  estimates.  Cox  used  this  approach  to  investigate  the  perfor¬ 
mance  of  the  HR  method  for  varying  signal  to  uncorreiated  noise,  and  to  discern  the  resolving 
power  between  two  signals  with  different  slowness  vectors.  The  HR  method  was  found  to  exhi¬ 
bit  higher  resolution  than  the  conventional  beamforming  method  in  the  precence  of  uncorre¬ 
iated  noise.  Interference  was  found  to  be  severe  when  steering  vectors  for  two  separate  signals 
are  perpendicular  to  each  other  in  the  metrics  defined  by  the  cross-spectral  matrix. 

Singular  value  decomposition 

Using  the  property  of  Hermitian  symmetry  for  S,  the  singular  value  decomposition  of  S 
and  S~'  can  be  written  as 

s -fxjV/V/*  (10) 

^■1 

s-*-f  xr'V/V,"  (id 

l-l 

where  (x,l  and  (V,]  are  the  sets  of  eigenvalues  and  eigenvectors  of  S.  Substitution  of  10  and 
1 1  into  3  and  9  yield  expressions  for  Pc,  and  PN* 

Pc-i>/U"V,V/<U  (12) 

P*m£x-,U"V(V("U1-1  (13) 

f-i 

The  two  methods  represent  different  weighting  schemes  for  the  eigenvector  contributions  to  the 
matrix  S  or  S~l.  The  formulas  in  equations  12  and  13  are  analogous  to  the  formulas  for  total 
resistance  of  a  network  of  resistors  in  series  or  parallel.  The  conventional  method  weights  the 
largest  eigenvalues  of  S,  while  the  HR  method  weights  the  smallest  eigenvalues  of  & 

In  the  case  of  a  single  plane  wave  of  amplitude,  nu\  propagating  across  the  array  with 
slowness  vector,  s'-k'/w,  the  cross-spectral  matrix  is 

J„-M  V,{\k)  Uj(k’)  (Id) 

where  rtf-(x,-xy,>/-y;).  The  matrix  is  singular  since  each  row  or  column  is  a  linear  combina¬ 
tion  of  any  other,  Often  in  practice,  S,  is  singular  or  nearly  so.  The  matrix  is 


prewhitened  with  magnitude  «,  and  inverted  with  standard  methods.  Prewhitening  as  suggested 
by  Capon  et  al  (1967)  can  be  accomplished  by  multiplication  of  off-diagonal  elements  by  1-c. 
This  is  equivalent  to  adding  uncorrelated  noise  to  each  station  and  makes  the  matrix  diagonally 
dominant.  This  is  equivalent  to  reducing  the  interstation  coherence  by  adding  uncorrelated 
noise  to  each  station.  The  prewhitening  serves  to  perturb  the  zero  eigenvalues  to  a  finite  value, 
so  X=*.  Symmetric  perturbation  of  the  Hermitian  symmetric  matrix.S,  produces  unpredictable 
rotations  of  the  degenerate  eigenvectors,  while  they  still  span  the  same  space  as  before  the  per¬ 
turbation  (Wilkinson,  196S;  Davis  and  Kahan,  1970). 

If  the  matrix  S  is  perturbed  to  the  matrix  S'  then  the  perturbations  of  the  eigenvalues 
may  be  described  by  an  expansion  in  the  basis  provided  by  the  old  eigenvectors, 

*  <  ▼>  *U1 

imjKi 

Obviously,  the  expansion  is  invalid  for  degenerate  eigenvalues  but  serves  to  explain  the  insta¬ 
bility  of  the  eigenvectors  for  pairs  of  eigenvalues  that  are  closely  spaced.  If  two  non-interfering 
waves  of  comparable  amplitude  are  combined  the  two  plane  waves  will  have  corresponding 
eigenvalues  of  comparable  size  and  render  the  two  eigenvectors  unstable.  This  confirms  the 
results  of  Der  and  Flinn  (1975)  that  signals  of  comparable  amplitude  may  be  difficult  to 
separate. 

In  the  case  of  a  sum  of  plane  waves  with  amplitudes,  and  slownesses,  s„,  the  cross- 
spectral  matrix  is 


ti,expliu  »„•!<,]+  £  l  05) 

m  n0m 

If  we  examine  the  decomposition  of  S  as  seen  in  equation  10,  the  decomposition  is  in  the  same 
form  except  for  the  cross-term  that  produces  the  interference  between  pairs  of  waves.  This 
observation  may  explain  the  difficulty  of  the  HR  estimate  to  resolv  wo  closely  interfering  sig¬ 
nals.  If  a  pair  of  signals  coincide  with  each  other’s  array  sidelobes,  the  estimate  may  be 
unstable  at  that  frequency  do  to  the  interference.  Smoothing  of  the  cross-spectrum,  SU>),  is 
optimal  when  the  interference  terms,  fluctuate  more  rapidly  than  the  simple  sums  of 


plane  waves. 


An  Example 

As  an  interesting  example  of  cross-spectral  decomposition,  an  array  experiment  conducted 
along  the  San  Andreas  fault  will  be  briefly  discussed.  Complete  details  of  the  experiment  are  to 
be  found  in  Chapter  4.  The  array  was  designed  to  measure  slowness,  azimuth  of  arrival,  and 
polarization  of  S  waves  from  earthquakes  within  the  fault  zone  3  km  away. 

Figure  17  of  Chapter  4  shows  a  high  resolution  and  a  conventional  beam  estimates  of  f-k 
spectra  at  6.6  and  7.8  Hz  for  the  east  component  of  motion  .5  sec  of  S  waves.  North  is  to  the 
top,  and  east  is  to  the  right.  Contours  are  in  decibels  with  respect  to  the  maximum.  The  two 
contour  plots  exhibit  some  of  the  common  traits  for  HR  and  beamforming  f-k  spectra.  The 
beamforming,  conventional  method,  shows  two  broad  maxima  coming  from  the  east  and 
southwest  quadrants.  The  HR  estimate  shows  much  higher  contrast  over  the  background, 
sharper  peaks,  and  a  more  complicated  background  20  db  down  from  the  maximum.  Superim¬ 
posed  on  each  plot  are  four  slowness  vectors  at  (.80  s/km,  110°),  (.55  s/km,  120°),  (.40  s/km, 
205°),  and  (.40  s/km,  250°).  The  four  vectors  correspond  to  four  beams  shown  in  Figure  17 
(Chapter  4)  for  horizontal  components  resolved  as  transverse  (T),  and  radial  (R)  along  the  four 
slowness  vectors.  The  four  beams  were  selected  from  a  set  of  beams  at  0.05  s/km  intervals  in 
the  general  areas  of  the  broad  maxima  of  the  f-k  diagrams  of  Figure  16  (Chapter  4).  The 
beams  represent  a  broader  bandwidth  and  therefor  are  better  estimates  of  the  location  of  impul¬ 
sive  arrivals.  The  beamforming  suggests  that  the  broad  conventional  f-k  maxima  are  multiple 
peaks  on  the  HR  spectral  estimate.  An  interesting  phenomenon  pointed  out  by  Cox  (1973) 
appears  to  occur  on  the  east  component  HR  f-k  spectra;  when  two  signals  of  different  strength 
arrive  with  closely  spaced  slowness  vectors,  it  is  possible  that  the  weaker  arrival  will  show  up  as 
stronger  on  the  HR  spectral  estimate.  This  is  clearly  a  disadvantage  of  the  HR  method  for  esti¬ 
mation  of  the  relative  strengths  of  interfering  waves.  The  time  domain  stacks  show  the  slower 
arrival,  .10  s/km,  1 10°,  to  be  stronger  on  the  transverse  component,  doeer  to  the  north  com¬ 
ponent  of  motion  for  an  arrival  from  the  east.  Conversely,  the  faster  arrival,  .55  s/km,  210°,  is 
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the  stronger  on  the  radial  component  and  should  show  nearly  equally  partitioned  between  east 
and  north  components  of  motion.  The  HR  f-k  spectra  show  the  two  arrivals  reversed  in  rela¬ 
tive  amplitude.  The  data  indicate  multipathing  is  important  in  this  geologic  environment  and 
anisotropic  S  wav  s$  may  be  responsible  for  a  significant  partitioning  of  the  two  arrivals  from  the 
direction  of  the  epicenter  to  the  east.  The  arrivals  from  the  southwest  may  be  birefringent  as 
well. 

The  data  serve  as  an  excellent  example  of  multipathing  and  interference  effects.  The  four 
arrivals  all  propagate  across  the  array  within  .5  second,  forming  a  complex  locus  of  travel  time 
planes  over  the  two-dimensional  array.  A  portion  of  the  multipathing  may  be  anisotropic  S- 
wave  propagation.  At  station  9  the  .80  s/km  arrival  is  first,  while  at  station  1,  the  .5S  s/km 
arrival  is  first. 

In  Figure  1,  the  eigenvalue  spectrum  is  shown  for  the  normalized  and  prewhitened  (a  fac¬ 
tor  of  .001)  cross-spectral  matrix  at  5.8  Hz.  The  cross-spectral  matrix  was  averaged  over  .8  Hz 
bandwidth.  The  eigenvalues  span  a  range  of  4  decades  and  without  normalization  and  prewhi¬ 
tening  would  span  15  decades.  The  cross-spectral  matrix  is  sorely  rank  deficient  with  a  condi¬ 
tion  number  of  10,000.  In  Figure  2,  the  HR  f-k  estimates  for  the  contributions  made  by  the 
smallest,  the  2  smallest,  the  3  smallest,  and  alt  the  eigenvalues  are  shown.  There  is  virtually  no 
difference  between  the  estimate  of  the  f-k  spectra  made  with  only  the  3  smallest  eigenvalues 
and  all  the  eigenvalues.  The  contribution  of  the  smallest  eigenvalue  is  clearly  the  most  impor¬ 
tant  and  corresponds  to  the  arrival  of  two  of  the  four  identified  arrivals.  Another  arrival  is 
given  by  the  next  largest  eigenvalue.  Interference  between  arrivals  produces  sidelobes  in  the 
slow  northwest  quadrant  and  to  the  south.  Comparison  with  the  conventional  f-k  estimate  of 
Figure  3  shows  the  diffuse  maxima  characteristic  of  the  beam  method  with  no  suggestion  of 
multiple  arrivals  from  the  east  or  southwest. 

Conclusions 

The  high  resolution  method  of  Capon  (1967)  differs  from  the  conventional  method  of 
frequency-wavenumber  estimation  by  utilizing  the  near  singular  components  of  the  cross- 
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spectra)  matrix  estimate.  Often,  the  cross-spectrum  contains  only  a  limited  amount  of  informa¬ 
tion  and  is  dominated  by  a  single  or  small  set  of  important  eigenvalue,  eigenvector  pairs  Pre¬ 
vious  analysis  of  the  HR  method  relied  on  either  a  multidimensional  noise  model  (Capon  1970) 
or  was  concerned  with  the  resolution  of  the  array  to  detect  and  separate  multiple  signals  of  lim¬ 
ited  slowness  separation  (Cox  1973).  The  present  analysis  shows  that  even  in  the  presence  of 
multipathing,  the  information  may  be  contained  in  a  limited  number  of  eigenvalues  and  each 
contributes  an  alternative  interpretation  of  the  phase  delay  data.  The  HR  method  weights  the 
inverse  contribution  of  the  near  singular  portion  of  the  cross-spectrum.  The  method  therefore 
gives  an  inverse  picture  of  the  null  space  of  the  data.  Consequently,  the  fictitiously  excited 
sidelobes,  and  sharp  peaks  of  the  HR  estimate  are  suggestive  of  real  signals,  but  must  be 
verified  by  additional  estimates  at  another  frequency  or  in  the  time  domain. 

A  method  that  would  yield  direct  estimates  of  the  maxima  for  f-k  spectra  without  detailed 
calculation  and  contouring  would  be  useful.  Such  a  method  would  require  finding  the  maxima 
of  the  inner  products  of  eigenvectors,  Vy,  with  the  station  phase  delay  vectors,  U(k).  Unfor¬ 
tunately,  it  is  evident  from  the  example  shown  that  a  single  eigenvector  may  contain  informa¬ 
tion  from  more  thar  one  arrival.  A  best  fit  to  plane  waves  may  not  be  the  only  interpretation 
for  each  eigenvector,  V y.  Furthermore,  interpretation  of  the  eigenvector  as  a  phase  delay  vec¬ 
tor  would  require  unwinding  phase  of  the  frequency  spectrum  in  two-dimensions.  A  sparsely 
occupied  two-dimensional  array  would  pose  a  formidable  computation  effort  for  uniquely 
unwinding  phase  in  two  dimensions  .unless  a  good  first  estimate  of  slowness  ts  provided. 

A  posability  suggested  by  the  analysis,  for  cross- spectral  matrices  that  exhibit  singular 
behavior,  is  to  solve  for  the  singular  eigenvectors  rather  than  to  prewhiten  the  matrix.  An  esti¬ 
mate  of  the  null  space  would  be  made  without  perturbing  the  matrix  and  its  eigenvectors.  The 
sharp  nature  of  the  inverted  null  space  is  not  yet  explained,  nor  the  fictitious  sidelobes.  Some 
of  the  results  of  Cox  (1973)  suggest  that  the  HR  method  may  show  only  gt eater  r.r>pare-.-*  ’■e-5 
lution.  The  location  of  the  sharp  peaks  are  not  more  precisely  located  than  the  sire  of  t! ;  t-vu 
lobe  of  the  beam  response.  The  HR  method  does  give  indications  of  closely  spw-<5  ,i  -  * 
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arrivals  not  apparent  by  conventional  means.  The  sidelobes  are  manifestations  of  the  interfer¬ 
ence  between  multiple  signals  that  can  not  be  unambiguously  removed  from  the  phase  delay 
data. 


81 

REFERENCES 

Burg  (1964)  Three-dimensional  filtering  with  an  array  of  seismometers,  Geophysics,  24,  S, 
693-713. 

Capon,  Greenfield,  and  Kolker  (1967)  Multidimensional  Maximum-likelihood  Processing  of  a 
Large  Aperture  Seismic  Array,  Proc.  IEEE  S3,  pg  192. 

Capon  (1970)  Probability  Distributions  For  Estimators  of  the  Frequency-Wavenumber  Spec¬ 
trum,  Proc.  IEEE  58,  10,  pg  1785. 

Cox  0973)  Resolving  Power  and  Sensitivity  To  Mismatch  Of  Optimal  Array  Processors,  JASA, 
54,  3,  pg  771. 

Davis,  G,  and  W.M.  Kahan  (1970)  The  rotation  of  eigenvectors  by  a  perturbation.  Ill  SIAM  J. 
Numer.  Anal.  7,  1,  pp  1-46. 

Der.Z.A.  and  Flinn,  E.A.  (1975)  The  applicability  of  Principal  Component  Analysis  to  Separa¬ 
tion  of  Multiple  Plane-Wave  Signals,  BSSA  65,  3,  pp  627-635. 

Jenkins  and  Watts  (1968)  Spectral  Analysis  and  its  applications,  Holden-Day,  San  Francisco, 
525  pp. 

Wilkinson  (1965)  The  Algebraic  Eigenvalue  Problem,  New  York;  Oxford  Univ  Press. 

Woods.  J  W.  and  Lintz,  P.R.  (1973)  Plane  Waves  at  Small  Arrays,  Geophysics  38,  6,  pp  1023- 


1041. 


•  f  * 
* 


82 


FIGURE  CAPTIONS 

Figure  1.  Eigenvalue  spectra  of  the  cross-spectral  matrix  for  the  east  component  of  motion  at 
5.8  Hz  with  a  0.8  Hz  bandwidth  for  the  short  S-wave  window  indicated  in  Figure  3 
(Chapter  4).  The  eigenvalues  span  4  decades. 

Figure  2.  Eigenvalue  decomposition  of  the  HR  f-k  spectral  estimate  of  the  east  component  of 
motion  for  0.S  sec  of  S  waves  at  5.8  Hz.  Ail  eigenvalues  were  used  to  estimate  the 
spectra  in  the  lower  right.  The  three  smallest  eigenvalues  were  used  to  construct  the 
HR  f-k  estimate  in  the  lower  left.  The  smallest  and  two  smallest  eigenvalues  were 
used  to  construct  the  two  upper  estimates.  Contours  are  shaded  to  -12  db  in  the 
upper  two  estimates,  and  shaded  to  -6  db  in  the  lower  two  estimates. 

Figure  3.  Conventional,  beamforming,  f-k  estimate  at  S.8  Hz  for  the  same  window  used  in  Fig¬ 
ure  8.  Contours  are  shaded  to  -6  db  w.r.t.  maximum. 
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IV 


THE  STOCHASTIC  GREEN’S  FUNCTION  AND  IT’S 
APPLICATION  TO  SEISMIC  SOURCE  INVERSION 


ABSTRACT 

Treatment  of  signal-generated  noise  (coda)  as  a  noise  source  for  deterministic  inversion 
of  seismic  data  can  be  formalized  by  the  uae  of  a  stochastic  Green's  function.  When  lateral 
variations  are  suitably  chaotic  the  random  fluctuations  of  elastic  waves  may  be  described  by  a 
stochastic  variation  upon  an  average  or  coherent  wave  field.  The  coherent  field  is  separated 
from  the  random  part  which  is  treated  as  signal-generated  noise.  Variance  estimates  for  the 
stochastic  part  of  the  Green's  function  can  then  be  used  to  estimate  source  parameter  variances. 
Signal-generated  noise  is  introduced  into  the  Green’s  function  deconvolution  for  seismic  source 
time  functions.  The  variance  estimate  of  the  moment  tensor  is  given  by  a  direct  application  of 
error  propagation.  An  example  is  given  for  an  explosion  source  at  NTS  using  estimated  seismic 
field  variances. 

Introduction 

The  simple  stratified,  layered,  approximation  to  earth  structure  allows  very  precise 
Green’s  functions  to  be  calculated.  These  approximations  often  contain  the  essential  wave  pro¬ 
pagation  characteristics  of  the  earth  structure  and  allow  inversion  for  the  seismic  source.  The 
lateral  variability  of  the  real  earth  makes  any  such  approximation  inaccurate.  In  a  situation 
where  the  earth  structure  is  dominated  by  the  vertical  variation  and  lateral  variations  are  suit¬ 
ably  chaotic,  a  laterally  averaged  structure  may  be  represented  by  a  traditional  horizontally 
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homogeneous  structure  with  a  superposed  horizontally  random  perM-y,  field  t**  r  * 
of  this  paper  is  to  demonstrate  that  the  problem  of  random  structure  is  separable  fro  n  'be  *>/e> 
agec  structure  wave  propagation  problem.  Green  function  estimates  *iveit  by  one-d:nw\  - 
structures  coupled  with  estimates  of  scattering  from  laterally  varying  structure  permit  the  treat¬ 
ment  of  the  scattered  signal  as  a  signal-generated  noise  process  for  t’.  e  of  source  inver¬ 

sion 

The  functional  formalism  for  the  solution  to  1iffer',r.t,iv  equations  vim  vmlom 
coefficients  is  reviewed  by  Adomain  (1964)  while  the  stv^vca1  "raiment  of  wave  propagation 
in  random  continuous  media  is  reviewed  by  Hoffman  *r„d  waves  scattered  by  oiscrete 

scatterers  is  treated  by  Twersky  (1964).  Keller  (1964)  reviews  some  of  the  approximations 
used  in  continuous  and  discrete  scattering  theory  and  demonstrates  that  under  suitable  condi¬ 
tions  transport  theory  (radiative  transfer)  is  an  equivalent  treatment  for  strong  multiple  scatter¬ 
ing.  Kara!  and  Keller  (1964),  Knopoff  and  Hudson  (1964),  and  Knopoff  and  Hudson  (1967) 
all  treat  the  vector  elastodynamic  equations  for  elastic  wave  propagation  in  a  continuous  random 
medium  with  special  attention  to  the  amplitude  and  phase  fluctuations  in  the  forward  direction. 
Frisch  (1968)  demonstrates  the  utility  of  diagram  methods  and  functional  integration  methods 
for  the  renormalization  of  wave  propagation  in  a  random  media  Dence  arc  Spence 
generalize  the  development  to  the  dyadic  stochastic  Green's  function  and  anisotropic  random 
media.  For  a  discussion  of  ensemble  averages  and  the  general  aprlir*tior  of  stochastic  random 
itdia  the  reader  is  referred  to  the  review  by  Hudson  (1982). 

The  effects  of  scattering  by  random  heterogeneities  includ-  amplitude  and  phase  fluctua¬ 
tions,  attenuation  of  the  average  or  'coherent"  held,  and  conversion  of  energy  betwet*  r=od?r 
of  propagation.  Formulation  of  a  dyadic  stochastic  Green's  function  has  the  advantage  of 
it  'j*  rg  all  of  these  affects  on  a  seismogram.  To  evaluate  the  inverse  problems  for  <*  ismu 
s  w  -cea  and  homogeneous  structures,  we  require  estimates  of  the  statistics  of  the  seism '  field 
i  ■  'n  the  random  variation  from  the  homogeneous  structure. 
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THEORY 

The  wave  propagation  problem  for  elastic  media  begins  with  the  linearized  equation  of 
motion  (Hudson,  1980), 


Lytij  —  ptii  —  (CyuUtJlj  “  /,.  (1) 

where  L0  is  the  elastodynamic  operator,  u,(r,t)  is  the  particle  displacement  at  the  position,  r, 
and  time,  r,  p(r)  is  the  density,  C^it)  is  the  linear  elastic  tensor  field,  and  /,(r,  t)  is  the  deter¬ 
ministic  source  field.  The  functional  dependence  of  the  elastic  moduli  and  density  on  position 
define  a  structure,  or  model  for  wave  propagation.  Perturbations  to  the  density,  p-p°+p',  and 
elastic  tensor  field,  C-C°+C',  are  introduced  such  that  the  elastodynamic  operator  may  be 
decomposed  into. 


(2) 

Lm  “  P°*w  “  (Cfjkfi **./)  j  (3) 

f-A  ”  P*®*  —  (CfjkAtk,l\j  (®) 

The  perturbations,  p1,  and  C1  are  considered  small  random  fluctuations,  <(p')J>  «  <p>2, 

<(C,U2>  «  <CJihl>\  and  <((CjL,)j)j>  «  <(CjUm)J>3,  and  assumed  to  have  zero 
mean,  <p'>— 0  ,  <Cy*fcl>-0,  <(CJlJ  J>-0.  The  operator  <  >  is  intended  to  symbolize 
the  ensemble  average.  Equation  (1)  defines  a  set  of  coupled  stochastic  differential  equations. 
Equation  (2)  defines  a  decomposition  of  the  stochastic  differential  operator  into  the  background 
operator  (equation  3),  and  the  random  part  of  the  operator  (equation  4).  Generally,  the  prob¬ 
lem  of  interest  to  seismology  would  be  the  case  where  the  average  density  and  elastic  constants 
are  functions  of  depth  only,  and  the  random  perturbations  are  functions  of  depth,  (z),  and 
lateral  coordinates,  (x,y), 

p(x,y,z)  -  p°(z)  +  p'(jr,.y,z)  and  C(x,/,z)  -  C°(z)  +  C'Oc, y,z). 

The  background  operator,  L £,  is  assumed  to  have  a  Green's  function  satisfying  the  equation, 

(5) 

with  the  appropriate  boundary  conditions,  such  that  an  inverse  operator,  (/£,  exists; 

LS,uf-f,  (6) 


implies  that 
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u,0(r,t)-G°/f  to 

-JVJV 1  *°(r,r\  t-fi 
The  solution  to  the  general  problem, 

Luuj  -  /,  (fc) 

“  //  — 

may  then  be  written  as 

«,  -  G°/y  -  G°Lj\uk  Vi) 

Expansion  of  equation  (9)  into  a  Neumann  series  (or  Bern  cxfans'  j-  1  produces, 

u,  -  (?«/,  -  G*L}kGlJ,  +  GlL\GlL}mGlJn  -  U0> 

—  uf  +  u }  +  u}  +  . 

Taking  the  expectation  of  equation  (10)  yields, 

<ut>  -  Gif,  +  GlKL^Ll^Glj*  +  .  HI) 

-  u,°  +  <  u,J>  +  . 

Conditional  on  their  convergence,  equations  (10)  and  (11)  serve  to  define  a  stochastic  Green's 
function  for  the  stochastic  differential  operator  Ltj  in  equation  (1).  Where  in  operator  notation, 

Gw  -  G°  -  Gl^Gl  +  GlL^GULG^  -  .  U2) 

<G„>  -  G%  +  +  (13) 

From  equations  (10),  (11),  (12),  and  (13)  we  decompose  the  displacement  field  two  ways, 

U,  -  U?  +  ur  -  Gif,  +  Gif,  04) 

ut  -  <«,>  +  ur  -  <G¥>f,  +  Gif,  OSi 

where  the  superscripts  *sc*  and  "st"  ate  intended  to  infer  ‘scattered*,  »r.d  *»ochastiu'.  I>  s  »;•  . 

from  equations  (14)  and  (15)  that  <■>  ^  a0,  and  <G>  G°  Note  that  the  exner.rr.e  • 

tally  determined,  or  averaged  field,  <u>,  does  not  represent  ii.s  background  structurr ,  -•  o 

and  <C>.  In  the  case  of  elastic  plane  waves  in  a  randomizer!  whole  space,  Karat  and  Keller 

(1964)  showed  that  the  presence  of  scattering  both  attenuates  and  disperses  the  expected  Vd 

rtlai.ve  to  the  background  solution,  u°.  The  measured  ve’oc-ue s  *:e  not  precisely  the  fere-- 

aioii  .d  structure’s  averaged  velocities  for  the  same  ttuc.i«  .h<  eKtxrted  field  if  so.t  »:*r 

solution  to  the  background  solution.  Therefore,  if  a  structure  is  experimentally  found  that 

gives  t  Green's  function,  G*«<  G>,  from  the  coefficients  p*.  and  CL\  then  it  does  not  f.-ibow 
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that  <p>-p*,  and  <C>— C*.  However,  this  is  just  what  is  usually  done  in  practice.  The 
laterally  homogeneous  structure  represented  by,  p*(z),  and  C*(i)  is  determined  that  matches 
the  observed  attenuation,  and  dispersive  effects  of  the  averaged  response,  <G>. 

We  will  decompose  the  Green's  function  into 

Gy  -  <  Gu>  +  Gfi  06) 

G"  is  the  stochastic  part  of  the  Green's  function  (Adomain;  1964;  Frisch,  1968;  Dence  and 

Spence,  1973).  Only  the  statistics  of  GJ'may  be  calculated  or  measured,  such  as  its  variance, 

<|G"IJ>.  Under  this  decomposition 

ut  «  <  G)j>fj  +  Gfjfj  m  <  U/>  +  uf  (17) 

and  the  stochastic  part  of  the  seismic  field,  u",  is  considered  signal-generated  noise.  The  vari¬ 
ance  of  the  stochastic  portion  of  the  seismic  field  is 

tw(u*)  «•  var( G"f) 

The  proportion  of  noise  in  the  seismic  field  may  be  measured  from  multiple  station  coherency 
estimates,  repeated  measurement  with  variable  source  receiver  paths  approximating  an  ensem¬ 
ble,  or  calculated  from  scattering  models. 

APPLICATION  TO  AN  INVERSION  PROCEDURE 

Stump  and  Johnson  (1976)  proposed  a  method  to  invert  for  the  moment  tensor  com¬ 
ponent  time  functions  of  an  equivalent  point  source  when  the  seismic  Green’s  function  is 
known.  For  layered  earth  models  precise  Green’s  function  estimates  may  be  calculated.  The 
presence  of  scattering  affects  an  inversion  in  three  ways  1.)  scattering  introduces  a  source  of 
unmodeled  attenuation  of  the  expected  or  coherent  waves,  2.)  the  signal-generated  noise  is  a 
fluctuation  of  amplitude  and  phase  of  the  expected  or  average  waves,  and  3.)  converted  waves 
may  appear  u  unmodeled  waves  producing  non-causal  model  source  terms. 

Attenuation  of  the  coherent  wave  field  by  the  generation  of  scattered  waves  has  recently 
been  approached  by  Aki  (1980,  1981),  Sato  (1981),  Kikuchi  (1981a,b),  Dainty  (1981),  Wu 
(1982a, b),  and  Sato  (1982).  Scattering  and  the  random  coda  noise  serve  to  introduce  uncer¬ 
tainty  into  the  amplitude  and  phase  of  the  observed  waves.  When  the  noisy  seismogram  is 
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deconvolved  with  the  appropriate  homogeneous  Green's  function  inttr-itc-i  .. 

source  time  function  estimate.  Amplitude  and  phase  fluctuations  of  the  expected  wave*  o. 
introduced  through  forward  elastic  wave  scattering  theories  such  a  Knopoff  and  Hudsu<.-  »  .< 
1967),  or  scalar  parabolic  theories  such  as  Chernov  (19o«l),  t  v-.i  ki  ‘,<977,1980),  or  McCoy 
(1980).  The  signal -generated  noise  due  to  coda  genercuon  -lay  «  rt.oduced  with  molc.s 
such  as  Wesley  (1965),  Dunkin  (1969),  Aki  and  Chouet  v'-97'.,  >  (1977  Va  •> 

(1978,  1980).  All  of  these  approaches  have  proposed  re;  ’  A  :■  r  :.c  calculator  of  th:  i  r/ 

of  the  stochastic  Green’s  function. 

Under  the  assumption  that  the  elastic  modulii  wsx  ,jc  .  >  are  random  function,  o i  *;,e 
horizontal  coordinates,  (x,y), 

<p(x,y,z)>  —  p°(z) 

<Cuk,(x,y,i)>  -  C,5i,.(z), 

the  background  problem  is  modeled  as  a  stratified  medium,  where  material  properties  vary  only 
in  the  vertical,  (z),  direction.  The  horizontal  heterogeneity  is  modeled  as  random  perturbations 
to  the  stratified  medium  and  the  decompositions  of  the  seismogram  and  Green's  function  into 
deterministic  plus  stochastic  parts  is  used  (equations  16,  and  1 ,). 

Following  Stump  and  Johns*.  <1976),  we  write  the  linearized  reiatior  for  d<  .V.; Place¬ 
ment  seismogram,  in  the  frequency  domain  as, 

U(/)-G(/)M(/}  N'/i-  1  1 

The  inverse  solution  for  the  source  becomes 

M’( /)  -  G-l( /'!’(/'  M9») 

wrtMVJl-G-'wrO'y.riS  v  9  « 

where  G  is  the  Green’s  function  matrix,  G-1  is  the  inverse  of  the  Green's  function  matrix,  and 

G~r  is  the  Hermitian  transpose  of  the  inverse  of  the  Green’s  function  matrix.  N(U  w  Uk 

ambient  noise  in  the  data  vector.  IK/)  is  a  vector  composed  of  complex  spectra  of  *>•<  win 

lowed  seismograms  at  the  frequency,  f.  And,  M(/)  is  nt  matrix  of  mot  pendent  exvecAi  u 

the  moment  tensor,  with  M’(/)  the  estimate  at  frequency,  i. 
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If  G— <G>+G“  in  the  presence  of  only  signal-generated  noise  (N(/)-4). 

U(/)  -  <G(/)>M(/)  +  G"(/'M(/)  -  <U (/)>  +  U*(/)  (20) 

If  we  approximate  <G>  by  G*  then  we  may  produce  an  estimate  for  the  moment  tensor  ele¬ 
ments  as 

M‘(/)  -  G*-'(/)U</)  (21a) 

-M(/)+G*-'(/)G”(/)M(/) 

vor<  M*>  -  G*-'wr(U)G*-r  (21b) 

-G*-'wr(G"M)G*-r 

The  seismic  source  estimate,  M’  is  contaminated  by  the  signal-generated  noise  term. 
Since  the  variance  of  the  stochastic  part  of  the  Green's  function  is  zero,  the  esti¬ 
mate  of  the  source  is  unbiased  in  the  frequency  domain,  and  it's  variance  is  related  to  the  vari¬ 
ance  of  G,  or  of  the  data,  U. 

AN  EXAMPLE 

Data  from  an  array  of  strong  motion  accelerometers  6  km  from  the  explosion  COL  WICK 
at  the  Nevada  Test  Site  will  be  used  to  illustrate  the  estimation  of  the  seismic  moment  tensor 
and  accompanying  variances.  The  signal-to-noise  ratio  is  estimated  directly  from  the  data  and 
used  in  the  inversion  for  the  seismic  source  spectra.  See  Chapter  two  for  a  description  of  the 
experiment  and  analysis  of  the  data. 

In  Figure  1  the  frequency-depe>.Jent  spectral  variance  var(U(/))-<  llK/)-U(/)l,> 
and  the  inferred  signal-to-noise  power  ratio  for  the  P  wave  and  P  coda  window  are  estimated 
from  three  stations  an  equal  distance  from  the  explosion.  The  three  stations  are  spaced  at  200 
meter  intervals  across  the  direction  of  propagation  of  the  P  wave.  The  acceleration  spectral 
amplitudes,  |U(f)l  for  the  three  components  are  shown  in  Figure  2.  The  signal-to-noise  power 
ratio  (SNPR)  was  estimated  from  the  inverse  of  the  normalized 
variance, SNPR~'  -  voriU(J))  *1  (//(/)  I"2.  Of  particular  note  is  that  the  variance  of  the  verti¬ 
cal  component  is  nearly  constant  from  1  to  10  Hz,  while  the  total  spectral  amplitude  declines. 
The  net  effect  is  for  the  SNPR  to  decline  from  5  to  10  Hz.  A  simitar  trend  is  observed  for  the 
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radial  component  spectra.  The  transverse  component  SNPR  estimate  fluctuates  around  1.  It 
was  previously  argued  (Chapter  one)  that  a  significant  portion  of  the  transverse  component  in 
this  time  window  prior  to  the  S-wave  arrival  was  scattered  from  the  P  wave.  We  can  conclude 
that  the  decrease  in  SNPR  is  a  result  of  the  proportional  increase  in  the  incoherent  signal 
energy  and  not  simply  do  to  a  decline  in  signal  strength. 

In  Figures  3A,B,C,  records  of  stations  1  and  4  are  compared  in  the  time  domain  as  broad¬ 
band,  and  bandpassed  signals  between  5  and  10  Hz  and  between  10  and  20  Hz.  The  two  sta¬ 
tions  are  only  200  meters  apart.  The  circled  number  is  the  normalized  maximum  positive 
correlation  between  the  two  traces,  The  different  response  of  the  two  stations  to  the  S 
wave  is  particularly  pronounced  (1.3  to  2.S  sec  after  the  P  wave).  The  mismatch  between  the 
two  stations  is  both  in  amplitude  and  phase.  The  correlation  values  estimate  the  signal-to-noise 
ratios  for  the  5-10  and  10-20  Hz  bands  as  1.2,  .83  for  the  vertical,  0.70,  0.40  for  the  radial,  and 
0.65,  0.50  for  the  transverse,  SNR'1  -  R^-l.  This  is  in  contrast  to  the  broadband  SNR  of 
2.5,  2.2  and  1.0  for  Z,R,  and  T.  These  signal-to-noise  ratios  reflect  the  average  over  a  5.0 
second  window.  The  initial  vertical  P-wave  motion  remains  coherent  in  the  10  to  20  Hz  band 
although  it  has  become  emergent  at  the  higher  frequencies. 

The  goal  of  this  section  is  to  illustrate  the  use  of  the  inferred  signal-to-noise  ratios  from 
the  COL  WICK  array  data  to  quantify  the  uncertainty  in  the  inverted  source  spectra.  Results  for 
the  vertical  component  P  and  P  coda  spectra  are  shown  in  Figure  4A,B.  The  acceleration  spec¬ 
tra,  the  twice  integrated  Green’s  function  spectra  for  the  explosion  source,  and  the  inferred 
explosion  source  are  plotted.  The  short  data  window,  1.28  seconds,  is  inadequate  to  control  the 
low  frequency  behavior  of  the  source  and  a  longer  window  enclosing  S  waves  is  inevitable.  A 
longer  window  also  incorporates  larger  quantities  of  "noise*. 

The  use  of  three  components,  and  a  longer  time  window  are  used  to  incorporate  the 
observed  S  waves.  A  S.12  second  time  window  wu  chosen  (Figure  3).  The  average  seismo¬ 
gram  and  its  variance  are  estimated  for  a  stack  of  seismograms  at  a  slowness  of  0.2  sec/km 
slower  than  the  P  wave.  This  aligns  the  principle  S  waves  on  the  three  components.  The  vari- 
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ance  for  each  component  is  estimated  from  comparison  of  the  complex  spectra  for  stations  at 
the  same  distance  from  the  event.  The  average  complex  spectra  are  then  estimated  as  a  linear 
combination  of  the  spectra  of  the  individual  stations.  The  Green’s  functions  for  the  different 
distances  are  similarly  stacked  at  the  same  slowness.  Therefore,  the  calculated  Green's  func¬ 
tions  have  been  given  the  same  spatial  filtering  as  the  data. 

Three  separate  estimates  for  the  isotropic,  or  explosive,  source  were  made  using  1.)  the 
radial  and  vertical  displacement  spectra  estimates,  (V2  and  U*),  2.)  the  radial  displacement 
data  only  (Ux),  and  3.)  the  vertical  displacement  data  only  (Uz).  Figure  S  shows  the  ampli¬ 
tude  spectra  of  the  Green's  functions  used.  An  overdetermined  explosive  source  estimate 
using  both  the  vertical  and  radial  component  data  is  shown  in  Figure  6.  The  square  root  of  the 
variance  (a)  is  plotted  as  a  dotted  line  in  Figure  6.  The  normalized  spectral  residuals  for  the 
radial  and  vertical  components  used  in  the  inversion  are  plotted  against  frequency  in  Figure  7. 
Since  only  S  seconds  of  data  was  deconvolved  by  5  seconds  of  Green’s  functions,  only  the  first 
S  seconds  of  the  resulting  source  estimate  may  be  considered  causal.  The  time-domain  far-field 
source  estimate,  in  Figure  8.  shows  non-causal  activity  between  S  and  10  seconds.  The  addi¬ 
tional  S  seconds  of  source  time  function  is  due  to  the  extension  of  the  data  and  Green's  func¬ 
tions  with  zeros  prior  to  the  calculation  of  the  first  discrete  Fourier  transform.  The  isotropic 
source  estimate  using  only  the  radial  component  data  is  presented  in  Figures  9  and  10.  The 
estimate  using  only  the  vertical  component  data  is  presented  in  Figures  11  and  12.  The  result 
based  on  the  radial  component  alone  is  the  least  non-causal  (Figure  10). 

Since  the  deviatoric  part  of  the  source  contributes  S  waves  to  the  vertical  and  radial  com¬ 
ponents  more  efficiently  than  an  explosive  source,  the  estimates  of  the  source  spectra  in  Fig¬ 
ures  6,  9,  and  1 1  are  conservative  upper  limits  for  the  explosive  source. 

SUMMARY 


We  have  shown  that  the  problem  of  elastic  wave  propagation  in  a  random  structure  is 
separable  from  the  propagation  problem  in  the  averaged  structure.  This  stochastic  propagation 
is  independent  of  the  source,  so  a  stochastic  Green’s  function  may  be  defined  (equations  12 
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and  13).  Under  the  decomposition  of  the  Green’s  function  into  it’s  deterministic  and  stochastic 
parts,  the  linear  inversion  of  seismograms  for  seismic  source  properties  remains  linear  and  the 
stochastic  response  of  the  medium  can  be  interpreted  as  signal-generated  noise.  This  noise  due 
to  scattering  may  be  empirically  estimated  or  derived  from  a  scattering  model.  An  example  of 
source  estimation  with  empirically  estimated  seismic  variances  was  presented  using  the 
COLW1CK  array  data. 
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FIGURE  CAPTIONS 

FI* art  1.  Smoothed  spectral  variance  estimates  a nd  Inferred  signal- to- noise  ratios  (SNR) 
for  the  vertical  (Z),  radial  (R),  and  tranaverae  (T)  component  spectra  of 
COLWICK  array  station*  1,2,  and  4.  1.28  second  window  encompassing  the  P 
and  P  coda.  The  0.79  Ha  bar  Indicates  the  smoothing  window  nsed.  The  loga¬ 
rithmic  scales  refer  to  the  SNR  estimates.  The  variances  have  been  multiplied  by 
an  arbitrary  scale  to  bring  them  onto  the  same  graph  as  the  SNR  estimates. 

Figure  2.  Acceleration  amplitude  spectra  for  the  1.28  second  window  used  in  Figure  1. 

Figure  3A,B,C.  5.12  second  broadband  and  bandpass  acceleration  records  for  the  COLWICK 
stations  1  and  4.  Broadband  records  are  on  the  left.  5-10  Hs  bandpass  records 
are  in  the  center.  10-20  Ha  bandpass  records  are  on  the  right.  The  value  In  a  cir¬ 
cle  is  the  maximum  normalized  cross-correlation  between  two  records.  Other 
numbers  above  and  below  are  peak  values  in  counts.  All  traces  for  a  given  com¬ 
ponent  (Z,R,  or  T)  are  the  same  scale. 

Figure  4.  A.)  (Above)  average  vertical  acceleration  amplitude  spectra  for  the  1.28  second 
window  from  COLWICK  stations  1,2  and  4  (solid  Une).  The  square  root  of  the 
modulus  of  the  variance  estimate  Is  plotted  as  a  dashed  line.  (Below)  Green’s 
function  (ramp  source)  amplitude  spectra  for  the  same  window  B.)  Isotropic 
source  (displacement)  amplitude  spectra  estimate  (solid  line)  and  square  root  of 
the  modulus  of  the  variance  estimate  (dashed  line). 

Figure  5.  Amplitude  spectra  of  the  vertical  (solid  line)  and  radial  (dashed  Une)  Green’s 
functions  for  the  Isotropic  (exlposlve)  source  appropriate  for  the  COLWICK  array. 
Green’s  functions  have  boon  stacked  on  the  S  wave  arrival  for  four  distances  4.00, 
4.087,  4.173,  and  4.344  km  from  the  source. 

Figure  4.  Estimate  of  the  far-leld  Isotropic  source  spectra  modulus  (solid  Une)  and  the 
uncertainty  estimate  (dashed  line)  derived  from  the  stacked  radial  and  vertical 
spectra  and  stacked  Green’s  function*. 
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Fit  arc  7.  NmmIM  mMnlt  of  tbo  vortical  and  radial  spoctral  components  a  sod  to  esti- 
mate  tbo  isotropic  uuct  of  Figaro  i. 

Figaro  8.  Tine  domain  plot  of  tbo  isotropic  source  estimate  of  tbo  source  estimate  of  Figaro 
d.  Considerable  non-causal  energy  arrives  between  5  and  10  seconds.  Units  are 

1020  <fyne- cm/  sec. 

Figure  9.  Estimate  of  tbe  far- filed  isotropic  source  spectra  modalas  (solid  line)  and  tbe 
uncertainty  estimate  (dashed  line)  derived  from  the  radial  data  only. 

Figure  10.  Time  domain  plot  of  the  far-ficld  source  function  of  Figure  9.  Units  are 

\0w<fyne-cm/  sec. 

Figure  11.  Estimate  of  tbe  far-field  isotropic  source  spectra  modulus  (solid  line)  and  the 
uncertainty  estimate  (dashed  line)  derived  front  the  vertical  data  only. 

Figure  12.  Time  domain  plot  of  the  far-held  source  function  of  Figure  11.  Units  are 

lOK<fyne~cm/  sec. 


LOGCDYNE-CM )  -20 


109 


L0GCDYNE-CM)-20. 


SCATTERING  OF  ELASTIC  WAVES  FROM  BOUNDED 


INHOMOGENEITIES  AS  A  MOMENT  EXPANSION 


ABSTRACT 

The  scattering  of  elastic  waves  by  an  inhomogeneity  in  a  uniform  matrix  is  addressed  with 
special  attention  to  incident  plane  waves.  The  solution  is  formulated  as  a  moment  expansion  of 
the  equivalent  force  and  stress  distribution  of  the  scatterer  due  to  the  incident  wave.  Although 
the  expansion  remains  exact,  rapid  convergence  is  only  assured  for  the  long  wavelength  regime, 
ka<l,  where  k  is  the  incident  wavenumber  and  a  is  the  scatterer  scale  size.  The  method  is 
shown  to  be  equivalent  to  the  form  function  or  Fourier  transform  approaches.  Conditions  for 
forward  conversion  of  P-to-S  waves  are  addressed.  Solutions  are  given  for  the  general  anisotro¬ 
pic,  Gaussian,  exponential  or  ellipsoidal  scatterers.  It  is  shown  that  trade-offs  exist  between 
models  for  the  shape,  heterogeneity,  and  anisotropy  of  scatterers,  but  that  a  general  scatterer 
may  not  be  modeled  by  a  homogeneous  isotropic  elastic  scatterer.  The  case  of  the  randomly 
heterogeneous  scatterer  is  treated  with  special  attention  to  the  random  Gaussian  scatterer. 
Under  the  conditions  of  ka<  <1,  application  to  the  attenuation  and  dispersion  of  the  effective 
wave  from  multiple  scattering  is  discussed. 

Feraralatioa  ef  the  PraMem-  the  Moment  Expansion 

Methods  are  required  to  describe  the  scattering  of  elastic  waves  by  small  scatterers  (small 
with  respect  to  a  wavelength).  Important  quantities  such  as  total  scattered  and  converted  elastic 
energy,  as  well  radiation  patterns  are  needed.  A  moment  expansion  of  the  scatterer  has  prom- 
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ise  for  the  calculation  and  analysis  of  scattered  wave  fields  from  small  bounded  heterogeneities. 
Several  methods  for  computation  of  scattering  cross-section,  and  radiation  pattern  of  scattered 
energy  for  compact  scatterers  in  elastic  whole  spaces  have  been  suggested  and  implemented  for 
incident  plane  waves  (Ying  and  Truell,  1956;  Knopoff,  1959;  Johnson  and  Truell,  1965;  Kraft 
and  Franzblau,  1971;  Mcbride  and  Kraft,  1972;  Pao  and  Mow,  1973;  Waterman,  1976;  Vasun- 
dara  and  Pao,  1976;  Gubernatis  et  al,  1977a,b,  Visscher,  1981;  Varadan  and  Varadan,  1982;  and 
Herman,  1982,  to  name  a  few).  The  methods  employ  either  special  function  expansions 
^spherical  harmonic,  or  elliptical)  of  the  incident  and  scattered  fields  or  a  solution  of  boundary 
integral  equations  of  the  Fredholm  type.  Common  draw-backs  of  these  methods  are  the  inabil¬ 
ity  to  treat  nonuniform  volumes,  irregular  shapes,  or  anisotropicaly  heterogeneous  bodies.  An 
equivalent  force  moment.  Green  function  derivative  expansion  or  moment  expansion  may  have 
some  advantages  as  a  description  for  general  compact  scatterers  if  the  object  lacks  symmetries 
necessary  for  special  function  expansions  and  if  the  multipole  expansion  converges  rapidly 
enough  in  the  frequency  bandwidth  of  interest.  The  use  of  equivalent  force  moments  to 
describe  an  indigenous  elastic  source  is  common  in  seismology  since  Gilbert  (1970)  introduced 
the  equivalent  stress  moment.  Backus  and  Mulcahy  (1976a, b)  explored  higher  order  moment 
tensor  expansions  for  indigenous  seismic  sources.  Higher  order  moments  for  propagating  faults 
have  recently  been  used  by  Stump  and  Johnson  (1982)  to  compute  synthetic  seismograms 
Backus  and  Mulcahy  (1976a)  point  out  that  the  moment  expansion  for  the  elastodynamic  source 

converges  at  the  rate  of  — r(-~)\  where  a  and  k  are  the  characteristic  source  size  and 

*!  k 

wavelength.  The  same  applies  for  a  moment  expansion  of  the  scattered  field  from  a  compact 
scatterer  of  size  a. 

The  general  equation  of  motion  for  the  elastic  body  can  be  written  as 

purlC&Htj),/*  ft  (1) 

By  either  introduction  of  Fourier  transforms  or  factorization  by  e*'  we  can  write  the  steady- 

state  harmonic  equation  as 


-Pw^HC^Wm),/-  ft 


(2) 


lid 


Let  the  inhomogeneity  be  represented  as  a  perturbation  to  a  problem  for  which  a  Oreen’s  func¬ 


tion  is  known.  The  density  and  elastic  variations  are  written  as 

p— p'+8p,  Cjy*/* C' yu+t  C^.  (3) 

The  background  solution,  u°,  is  the  solution  to 

ft  (4) 

and  the  perturbation  or  scattered  part,  u*  is  the  solution  to 

-p^uf-iC'uuulJj  ~  8 pwJ«,  +  ( SC0UukJ)  j  (5) 

where  u ,  is  the  complete  solution,  u,  -  u,°  +  uf.  If  the  Green's  function  exists  for  our  back¬ 
ground  equation  (4)  then, 

fjg0(t,t')  </V  (6) 

where  R°  is  the  region  bounding  the  equivalent  body  forces,  fy  We  then  have  that  the  scat¬ 
tered  field  is  given  by 

n5»)“,/jr(*P«Ju/r')  +  (8  CJklm  (f,)uim(t'))%k]gu(  r,  r’)  dV  (7) 

where  R’  is  the  region  bounding  the  inhomogeneity.  Using  the  divergence  theorem  and 

integration  by  parts, 

uflr)  -  JiJIl8pwJu;(i/)fv(r,rO  -  8  Cj^  ( rO  u, r.iOldV  (8) 


The  latter  expression  contains  no  derivatives  of  the  inhomogeneous  elastic  constants  and  only 
first  derivatives  of  the  total  displacement  field.  Less  must  be  known  about  the  inhomogeneities 
to  specify  the  scattered  field  than  in  equation  (S).  An  equivalent  view  can  be  derived  if  we 
define  the  equivalent  scattering  force  distribution,  flit), 

fit)  -  8p(r)«su,(r)  +  (8C0*,(r)u*,/(r))v/  (9) 

and  the  equivalent  stress  distributions  are  those  fields,  <r£(r),  such  that  /'-  -<r  The 

nonuniqueness  of  such  equivalent  stress  distributions  is  well  known  (  Backus  aod  Mulcahy; 

1976a,b).  The  formalism  remains  valid  for  inbomogeneities  incorporating  discontinuities  of 

displacement  or  stress  such  as  cracks,  voids,  or  fluid  inclusions.  <ry(r)  is  replaced  by  the 

equivalent  stress  field  that  would  produce  a  stress-strain  field  upon  a  surface  surrounding  the 

inhomogeneity.  An  example  of  such  a  procedure  is  Eshelby’s  (1957)  static  solution  for  eUip- 
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soidal  inclusions  by  an  equivalent  stress  tensor.  Gubernatis  and  Domany  (1979)  use  Eshelby’s 
results  as  a  quasistatic  approximation  for  the  equivalent  stresses  produced  by  long  wavelength  P 
and  S  waves  incident  upon  ellipsoidal  flaws.  The  equivalence  relations  follow  from  the 
representation  theorem  of  Burridge  and  Knopoff  (1964). 

Ji|f,*(*.*')/*,(*')rfV  - 

JM.lg'*(*x')<r<»(*')»/  -  - 

Jtn  8*  <*.*')  <T  fad2*' 

g  is  the  Green’s  function  in  the  absence  of  the  inhomogeneity,  g’  is  the  Green’s  function  in  the 
precence  of  the  inhomogeneity,  fi  and  <r  L  ere  the  scattering  equivalent  force  and  stress  distri¬ 
butions  on  the  scatterer  surface,  9R',  and  <r£  is  the  equivalent  stress  distribution  on  the  sur¬ 
face,  9R ,  surrounding  the  scatterer. 

It  has  become  customary  in  recent  years  to  expand  solutions  in  terms  of  force  moments 
and  Green's  function  derivatives  in  the  form  of  a  Taylor  expansion  of  equation  (7). 


uHr)  -  k(()g0,kl  *,<{,r) 

«-o 


(10) 


where  the  force  moments  are  defined 

yi?  *,„,<<>  -JJ,/,<r')(rVf*1) . (11) 

for  a  point,  {,  the  centroid  of  the  equivalent  scattering  force  distribution,  ff, r').  In  the  case  of 
indigenous  sources  with  no  net  torque  and  no  net  momentum,  yj0)-0,  and  yjt’-yj)’.  How¬ 
ever,  the  'pseudo*  forces  that  represent  the  scattering  are  only  part  of  the  tout  solution  and  do 
not  necessarily  comply  with  these  conditions.  The  zero’th  order  moment  may  not  be  zero  and 
the  Arst  order  moment  may  not  be  symmetric.  Now,  it  should  be  easier  to  deal  with  the 
equivalent  stress  than  with  equivalent  stress  derivatives  so  by  use  of  Green’s  theorem  once 
again  we  may  write  the  moments  so  ss  to  eliminate  the  derivatives  of 

rl0)  ■  JM,  8pwJ«,dJr  +  JKf  (8  CyuUkJjfr  (12) 

-  Jj,  9pw'u,fr  +  Jm(  tCuuUtjHjfr 

where  i *y  is  the  normil  vector  to  the  surface  9R1.  If  the  volume  R'  is  chosen  large  enough 


such  that  the  integrands  go  to  zero  on  the  boundary,  then  the  zeroth  order  moment  simplifies 
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to 

y}0>  •  (13) 

We  have  chosen  a  coordinate  system  such  that  f  •  0.  Similarly,  for  n>l, 

yi"’...*.  -  *2fM,  tPK'k,-  (H) 

“  SK,  .  r^r. 

The  moments  depend  only  upon  the  inhomofeneity  fields,  gp,  &CIJkh  the  displacement  field,  u„ 
and  the  strain  field,  ui  js  within  the  scattering  volume,  R’.  The  moment  expansion  is  exact  as 
stated,  but  approximations  will  surely  be  required  for  application. 


The  Fourier  Transform  Approach 

We  may  compare  the  solution  given  by  equations  (10)  and  (11)  with  another  formalism 
used  in  the  fields  of  electromagnetic  and  quantum  scattering.  This  formalism  will  be  referred 
to  as  the  Fourier  transform  approach,  and  requires  the  use  of  the  far-field  approximation.  We 
will  let  our  Green's  function  be  that  for  an  isotropic  homogeneous  elastic  whole  space  for 
Irl » Ir’l  so  that  we  may  write 


f</(r,r')  88  g^tr.r')  +  f*(r,r’)  (15) 

,?'y)  • 

where  r  is  the  unit  vector  in  the  direction  of  r,  and  ko,  k0  are  the  P  and  S  wavenumber  vectors. 
The  far  field  scattered  wave  is 


Ignoring  all  terms  of  order  (~)2»  or  higher,  in  gv>k  equation  (16)  becomes 

uflr)  -  /;(k«,u)g^r)  + 
where  k,  -  kmr,  kt  -  kfr,  and 

fj(ki,u)  -  -  VcrkiCjHl,u,,m)*J9(-lk(r*r'))dir' 

We  recognize  /y(k)  as  the  3-D  spatial  Fourier  transform  of  the  equivalent  scattering  force 


(17) 


(18) 
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distribution.  The  results  of  equations  (17)  and  (It)  tre  equivalent  to  the  results  of  Gubernatis 
et  al  (1977),  generalized  from  a  homogeneous  scattering  body  to  a  general  body. 

In  order  to  see  the  connection  of  equation  (17)  to  the  moment  expansion  of  equation 
(10)  we  expand  <up(ik*r)  in  a  Laurent  expansion  within  the  integral  of  equation  (It),  k  —  *t. 


exp(iktr',)  -  £  —  1 


(19) 

(20) 


//« "  •  •  •  r'*A' 

We  recognize  the  integral  in  equation  (20)  as  the  moment,  >  j*** . . .  *<(  and  note  that  in  the  far 
field 


(r<f')  (-*J*?a,  •  •  firS)  (21) 

where  represents  differentiation  at  the  source  point,  r'.  The  far  field  Green’s  function  is 
symmetric;  Consequently  we  can  write 

«,*(')  *  7y(k«)^r)  4-  /y(k^)^r)  (22) 

*  Id/"!)  r ;•>..*,(  ...*•>)  +  4*-  *.(r)J 

M.^r)  -  £(l/"!>  k(r). 

The  first  two  expressions  are  approximate;  they  contain  only  far  field  terms  while  the  third 
expression  could  contain  all  terms  if  desired.  Furthermore  it  follows  that  the  first  moment 
terms  (n**0,l)  are  equivalent  to  the  long  wavelength  approximation  Gca<<l,  a  the  scatterer 
scale  length)  where 

7y(k)  =*  «yj  Bp  rfV  4-  <*?»«,.  J  8C**  rfV  (23) 

,  cxp(-/k*r')=l  ,  and  Uj  and  ity,*  are  nearly  constant  over  the  integration  volume  for 

k*r'«l.  The  Rayleigh  (long  wavelength)  and  far  field  approximations  then  yield 

u/*  uJmif&pdir'gIJ(r,r')  (25) 

-  tkUimf tCju.fr'lgljk.  +  g$kj 

-  -  Uj*2f  Bp rk'd3r'(g'km  +  fjk*). 

This  form  for  the  scattered  field  has  been  used  by  Knopoff  (1959),  Miles  (I960),  and  Guber¬ 
natis  et  al  (1977b,c),  usually  in  conjunction  with  the  Born  approximation,  «  *  u*.  The  last 
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term  of  equation  (25)  has  often  been  ignored  since  it  contains  an  additional  power  of  frequency 
over  the  first  term  involving  ip. 

The  first  degree  moments  are  long  wavelength  limits  of  scattering  strength.  The  higher 
moments  contain  information  about  the  orientation  of  the  scatterer  and  its  shape  and  sym¬ 
metries.  Care  should  be  exercised,  noting  that  the  nth  term  of  the  expansion  is  proportional  to 
(fcfl)" 

— — — ,  where  a  is  the  characteristic  size  of  the  scatterer.  The  expansion  may  converge  slowly, 
or  not  at  all  for  high  frequencies.  In  the  long  wavelength  limit  the  moments  (n>l)  become 
y i"*.  *.  =  w1",/ •  •  •  r’*,dV  (26) 

“  ul,mf  &Cikxlmr’ '  '  '  r>  km^r'' 

The  Forward  Scattering  Theorem 

Gubernatis  et  al  (1977a)  showed  that  the  optical  forward  scattering  theorem  may  be 
extended  to  an  elastic  body  in  a  whole  space  if  the  interference  between  converted  waves  is 
accounted  for.  The  scalar  forward  scattering  theorem  states  that  the  toul  cross-section  of  the 
scatterer  (integrated  over  all  scattering  directions)  is  proportional  to  the  phase  delay  of  the  toul 
field  observed  in  the  forward  direction.  This  theorem  is  of  great  importance  because  the  for¬ 
ward  direction  phase  delay  is  often  easier  to  measure  or  compute  than  the  toul  cross-section 
integrated  over  all  directions.  The  case  of  elastic  waves  in  an  isotropic  whole  space  requires  an 
additional  term  from  the  interference  of  scattered  P  waves  from  incident  S  waves.  The  far-field 
scattered  P  and  S  waves  are  given  by  equations  (22)  and  (21)  as 

,  A,{i)expUkmr)  £,(r)ejq>(tfc«r) 

If/  —  -  +  - c — ,  (27*) 

r  r 

In  terms  of  moment  expansions  A(?)  and  B(t)  are  given  by; 

A,(\)  -  — j-  I(l/i.!)(^«rtl--- (27b) 

«irpa  „ 

*,(«  -  — ^  I(l/n!)(V’?a,  •  1 '  U7c) 

The  incident  P  and  S  plane  waves  are  of  the  form 


if,0-  a8|,e*p(&>X|)  +  b,exp(lk0Xt)  0. 
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We  define  -i/(d). 


A,C)  - 

where  9  is  the  polar  angle  with  the  xt  axis  and  4  >*  the  azimuthal  angle  in  the  x*  xj  plane. 
The  theorem  as  stated  by  Gubernatis  et  aL  (1971a)  is  that  the  total  normalized  scattered 
power,  <r,(w)  is 

<rltt  -  / d(l  r^iCyuuljuj  (28a) 

where  / •  is  the  incident  power,  J dfi  is  the  integral  over  the  total  solid  angle  and  or,  is  the 

tout  cross  section. 


. .  ,  .  [«^i(0)«‘  +  «*,**(*)*;+ /5JA*(0)b; 

■,(»)  -  -  4»/««^ - JiiF+Tili* - 


(28b) 


cos(fl)  -  kjkf.  The  29  cone  describes  the  surface  where  interference  between  the  scattered  P 
wave  and  the  incident  S  wave  is  constant.  A  similar  cone  for  scattered  S  waves  interfering  with 
incident  P  waves  does  not  exist  so  long  as  a>0.  Anisotropic  media  with  psuedo  P  and  split  S 
waves  would  have  additional  terms  corresponding  to  similar  interference  cones. 


If  we  examine  equation  (28a)  for  an  isotropic  elastic  whole  space  then  upon  substitution 
of  equations  (27a,b,c)  we  find  the  total  scattered  power  may  be  given  as  an  expansion  of 
moment  moduli  and  integrals  over  the  radiation  patterns. 

a,r  -  / rfO—wplAl*  +  (2tc) 

If  we  examine  the  expansion  of  the  integrals  with  equations  (27b,c)  then  we  will  encounter 
terms  like 


Xj  . 

I 

zf  . 

I 

For  some  decompositions  of  y/0)  and  y™  the  integrals  are  given  in  Appendix  A. 


Forward  ceaverslea  ef  P-te-S  seder  the  Bern  eppresleutlen 

Under  what  conditions  can  we  obtain  P-to-S  conversion  in  the  forward  direction?  If  the 


incident  plane  P  wave  is  given  by. 


u?-b,\«xp(ikmx  |) 

then  we  must  get  a  non-zero  radiation  pattern  for  some  term  of  the  expression  for  the  scattered 
S  wave; 

utr)  -  yj0)(«  <,-?/,) 5(f)  +  y<(/urk(6irrtrJH-ikJS(r) 

+ . 

The  first  term  requires  for  j— 2,  or  3  that  y  j*y\.  The  requirement  for  the  second  term  to  be 
non-zero  is  that  the  principle  axes  of  y  do  not  coincide  with  the  x(  axis.  Examination  of  the 
form  of  y  ,(0), 

y<0)  -  ^ftpu^r 

,  makes  it  clear  that  without  strong  scattering  yj»  yj-  0.  If  the  scatterer  is  nonsymmetric, 
even  under  the  Born  approximation,  the  first  order  moment  may  be  of  the  form, 

y^'W/gpu^A  -  JbC^u^r 

-  (a,bj  +  OjbJ+Cj 

where  a  or  h  (a*b-0)  have  a  component  in  the  x  direction.  Consequently,  under  the  Born 
approximation,  the  zeroth  order  density  perturbation  does  not  produce  P-to-S  conversion  while 
the  first  order  elastic  perturbations  may.  If  the  elastic  constants  are  isotropic,  then 

(8k  d</+2  btffr 

The  first  order  moment  contributes  forward  P-to-S  converted  waves  only  if  the  scatterer  is 
nonsymmetric  about  the  x  axis.  Certainly  if  the  elastic  constants  are  anisotropic,  then  even  a 
symmetrically  shaped  scatterer  may  convert  P-to-S  waves  in  the  forward  direction. 

The  Gaussian  and  Exponential  Scattmri 

Let  the  density  and  elastic  constant  perturbations  be  proportional  to  an  exponential  distri¬ 
bution, 

bpexp(-r/a),  8  C^sjpf-r/e). 

Then  the  first  few  moments  under  the  long  wavelength  and  Born  approximations  are 
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yl0)-tw  aWbp  (29) 

ri‘)-8ira,«,t(I,dC**, 
y  Ju  "32  w  a* w3buufip 
y$L,m32nasbk,uf%'bCliim 
Or  if  the  scatterer  is  assumed  to  have  a  Gaussian  distribution 

bpexpi-rVai*)  bC^expi-rVa3). 

Again,  given  the  long  wavelength  and  Born  approximations,  the  the  first  few  moments  are  of 
the  form... 


<4)  -  *3  .,7  1/2  I 


y  tklmrt 


16 


y  i0>~a}ir  i,3u3ufip 
yV-ain>'\m6ClUm 

y  tu»aiY*in*1*uufip 

y  iuL-o5  |V/J#  to, C*, 
if 


(30) 


-p  if  (k>-/^m—  n,  or  ,Ar— or  /, 


otherwise,  zero. 


y%L,-  j|-eV/J|  some  05  y(4)] 


The  Ellipsoidal  Scatterer 

Ellipsoids  have  often  been  used  as  models  for  elongated  inhomogeneities.  The  long 
seismic  wavelength  compared  to  such  geologic  lenslike  bodies  makes  the  Rayleigh  approxima¬ 
tion  particularly  appropriate.  The  moments  for  a  homogeneous  ellipsoidal  scatterer  with  axes 
(01,0*0))  in  the  x  y  i  directions  has  moments  of  the  forms 

yj^-y-irfliOjflyu^yip  (31) 

yjkl-  -jy#  0,0,0  J mHb  kiO.a-b  *8  *a,o,)  ufip . 
yjiL—  ^y*o,e,  o,(«  *.0,0,-#  ^8^0,0,)  ii,>#8  Cm. 

Since  the  higher  order  terms  decay  so  rapidly  for  ka<<l,  density  contrasts  have  little 
affect  on  any  directional  properties  of  the  scatterer  for  long  wavelengths.  The  equivalent  stress 
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due  to  elastic  moduli  contrasts  are  the  dominant  source  for  any  anisotropic  scattering  For 
example,  consider  a  homogeneous  sphere  with  a  difference  in  shear  moduli  between  the  host 
medium  for  shear  in  the  xj  direction;  8C|j|j-8Cij2j-l/28Cj3jj-«.  If  no  density  contrast 
exists  then  the  first  moment  of  significance  for  an  incident  P  wave  with  wavenumber  vector  k 
is, 

yjJ>--jtfK«(8„  +  8*J) 

P  waves  traveling  in  the  xi  or  x2  directions  are  unaffected  by  the  scatterer.  Without  an  aniso¬ 
tropic  set  of  moduli,  a  homogeneous  inclusion  can  not  be  responsible  for  Rayleigh  scattering. 

Trade-Offs  Between  Shape,  Inbomogenelty,  and  Anisotropy 

Although  the  ellipsoids  and  Gaussian  scatterers  are  convenient  models  of  voids,  cracks, 
and  defects  a  strongly  inhomogeneous  scatterer  such  as  a  partially  filled  void  has  some  interest¬ 
ing  properties.  Consider  the  case  of  a  sphere  with  two  hemispheres  of  different  properties, 

8p'  tCjkh,  ,z<0,  and 

ip2  BCjitm  ,*>0. 

gp-fgp'+gptyz,  Ap-fsp'-sp7),  ac-tec’+ac1)/^,  Ac-lac'-ac1) 

The  differential  motion  of  the  top  and  bottom  of  the  sphere  even  under  the  long  wavelength 
and  Born  approximations  will  yield  equivalent  stresses  that  produce  scattered  S  waves  in  the 
forward  direction  from  incident  P  waves. 

raWujip  (32) 

yJi,^wa4ti2(A8p)8)tu/  - 
yffl-jj-v  aWbutijip 

+  J-(8ilg„+8ri8,2) 

+  yfi»l®ijlA8p  4  -90*ti,ufn^bCm 

ypL— 

+  *8  «*(-y-8  *8,1+808,1)  4-  j8,j8  J u,  fA8C^, 

+  ^■«‘w,8,*8  *8^(8  n8  *1+8  08^28^8  ,ol  «yA8p 
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Clearly  a  trade-off  exists  between  shape  of  a  scatterer,  an  inhomogeneous  seatterer,  and 
an  anisotropic  scatterer.  For  example,  the  partially  filled  sphere  has  moments  resembling  those 
of  a  homogeneous  anisotropic  ellipsoid.  Similarly,  a  homogeneous  sphere  filled  with  anisotropic 
material  may  resemble  a  homogeneous  anisotropic  ellipsoid.  However,  an  anisotropic  or  inho¬ 
mogeneous  scatterer  can  not  be  generally  modeled  by  a  homogeneous  scatterer  of  arbitrary 
shape.  For  example,  a  penny  shaped  crack,  void,  or  fluid  inclusion  exhibits  anisotropic 
behavior  and  can  not  be  modeled  by  a  homogeneous  scatterer  of  any  shape. 

The  Randomly  Inhomogeneous  Scatterer 

Suppose  the  density  and  elastic  perturbations  within  the  scattering  volume  can  be  con¬ 
sidered  random.  The  random  fields  are  defined  by  the  statistical  moments  where  <8p>-0, 
<80-0,  and  /tp„-<8p(r)8p(r')>,  are  their  spa¬ 

tial  auto-  and  cross-correlations  of  the  perturbed  density  and  elastic  fields.  The  expected 
moments  are  given  by, 

<y  ,<0,>— <8 put>  d3r  (33) 

<>**’  *(f>-wJ/<8pu/rft)...r4>>d,r 
—  J  <&Cn,iimui,mrki '  '  ' 

The  expected  scattered  field,  <  uf>,  is  given  by  an  expansion  of  equation  (10)  with  the 
moments  replaced  with  their  expectations.  The  moments  for  a  centered  random  field  will  be 
zero  in  the  Born  approximation.  The  variance  of  the  scattered  field  for  < «*>— 0  is  given  by 

m,  m 

If  we  ignore  the  cross  component  terms  ii*j)  then  we  have 

<U/l,>-Xl/(«!)<lrpi,...*.l,>!lw.*i  O*4"  moment  mrm  (34) 

« 

The  cross  moment  terms  may  be  significant,  but  we  will  presently  examine  only  the  main  corre¬ 
lation  terms. 

< ly/W,l*>— «*4<  / •pu/d’rj  fs pii,'rf,r2>  no£ 


(35) 


<  ly i°l2>-w4< J lf>utTkfrx  f Bpu'r^rj  no£ 

* 

+  2 Re(<  fipu.r^rt  J bC^ul^r^  >) 

If  (he  displacement  and  strain  fields  are  uncorrelated  with  the  random  inhomogeneity  fields, 

<ly/W)|J>*w4/  rfV|ii,|lJ  *„(r)d»r  (36) 

<  lyi,>tJ>— w4J*  J t0p(r)rkd*r 

+  Rcc-.J')S, 

+  2Re(/  R'C^fr) 

In  the  long  wavelength  approximation 

<ly/0,l2>s=«4l«<lJJ*  d*r'f  (37) 

<  ly*l)lI>  ^u*\uyj  (fir'f  /?P),(r,r')r*dJr 
+  “/,  *.«,'«/  fr'f  *cc'lk^r'r')dSr 
+  2Re(ulu,.tJ  dVj*  RpClUm^y)rk<fr) 

The  relationship  between  the  moments  of  a  distribution  and  its  derivatives  of  its  Fourier 
transform  can  be  used  to  facilitate  the  computations.  Following  Bracewell  (1978)  we  write 

yjkl  *.  ’  •  •  rk/j(r)^r  (38) 

(— 2w/)" 

where  t<(0)  is  the  nth  order  partial  derivative  of  the  3-D  Fourier  transform  evaluated  at 
k-0.  Since  the  auto*  and  cross-conelation  are  Fourier  transforms  of  the  auto*  and  cross¬ 
spectra,  the  nth  order  moments  may  be  quickly  interpreted. 


l8p(0)|J 


J*p,(r)i>d»r 
J*„(r)  r;Vr- 


1  d8o(0) 

2w/  6kj 

i  a28o(o) 

(2pt)*  bk} 


In  the  case  of  a  Ouassian  random  field  for  density  and  elastic  moduli  with  characteristic 


lengths  01,02,0), 


fi(f)  »  txpl-w(rf/ai+ r\ta\ +r}/*})l 
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h(k)  —  expl-w(a}kf+a}kt+a}k})) 

R„„(r)  -  <8p2>/i(r),  R^.- <btbt>  h(r),  RfC- <bpbC> h(r) 

K0-4*ai«2«^3 

And  the  flrst  order  terms  for  scattering  volumes,  V>  >  V0  ,  *re 

<l>?ll>=s  «/|u<lJK<8p,> 

<!yi,,l2>=s  w4lu,l2K<8pJ>e* 

+  u,%mu','tV<bCbC>  tu^  (39) 

+  2  Re(u,uimV<bpbC> 

If  the  inhomogeneity  field  is  statistically  isotropic  with  Gaussian  correlation  length,  a,  statisti¬ 
cally  isotropic  real  elastic  moduli  leads  to 

<lyj*,i2>a=  w4lu;l2«*2F<8pI>  +  |uu|JK<8A,>8>k  (40) 

+  +  4|u>t/|J  K<8X8p>8j* 

Unless  <8p8C>  is  imaginary  the  cross-correlation  between  density  and  elastic  moduli  does  not 

contribute  to  the  first  two  terms.  We  can  write  down  the  expected  energy  radiated  from  the 

random  Gaussian  scatterer  as, 

<  P>  —  Imag(.r2J d(l  fjC uUultlUj')  (41) 

+  ^-ap(— ! ~)>lc^(/^+/^)+/*G) 

2  4irpa 

♦  ^M^W.UtE+hD+ltG). 

£<|r»)|  »>-/• 
j 

<!y»,|J>-£8(/+J>cj+Gdj 

where  Dmc^  it  the  equivalent  dislocation  and  Gx,2dv  is  the  remaining  deviatoric  pert  of  the 
scatterer. 


For  a  P  wave  incident  with  amplitude,  bh  and  wave  number  vector,  kh  the  tout  scattered 
power  may  be  computed  from  equations  (41)  and 

f-«4F<  gpW  (42) 

E  -  lkl,»'(<•A,>♦2<•p,>44<B^a*>♦2<•p,>)l^, 
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0-0 

Gdjit  —  **u}a}V  <8p*> . 

For  a  plane  S  wave  incident  with  amplitude  vector,  v„  and  wavenumber  vector  kh  the  total 
scattered  power  may  be  computed  from  equation  (41)  and 

F-*4K  <«pI>vJ 
£-0 

l*yv*+**VylJK<aMJ> 

Gdjl  —  w4v/a*  F<6pJ> 

To  obtain  the  cross  section  and  hence  the  attenuation  of  the  incident  wave  we  require  the 
normalization  of  P  by  the  input  power,  /*, 

/•-  ^lapluP+ZipIvi1].  (45) 

For  an  incident  P  wave  <af>»<P>/(a>iap),  and  for  an  incident  S  wave 
<<t ?>~<S>/ (w2/Jp) .  From  Twersky’s  (1964)  multiple  scattering  scalar  wave  theory  for 
point  scatterers,  the  attenuation  of  the  effective  wave  that  propagates  through  a  random  distri¬ 
bution  of  the  point  scatterers  with  density  p,  is  given  by,  *o(l+v),  P—2 irQ/ki,  and 

Imag(v)-p  Jl.  From  equations  (41),  (42),  and  (43)  we  can  calculate  the  attenuation  of  the 
effective  P  and  S  wave,  lmag(Kt//),  for  a  random  distribution  of  random  scatterers.  The 
effective  P  and  S  wave  turbidity  coefficients  T*+rK,  and  rc+T^  are  given  by 

r'V  -  p,(— ^T)2l/,£+(/2£+/^)+/«<7)*iJ 

4irpo* 

r*V  -  p,|-(~^3-),[/s£+(/*£+/70+/«G)*jI 
and 

TvvJ  -  p.^—^WliF+UiE+hD+ItOki] 

P  4wpaJ 

-  p Sj^jpWhF+UtE+hD+ltOk}) 

with  F.E.D,  and  O  given  by  equations  (42)  and  (43)  respectively.  Care  must  be  exercized  in 
evaluation  of  Re  (£,//)  since  the  forward  scattering  intensity  9  contains  a  term  due  to  the 
interference  of  P  and  S  waves  not  evaluated  in  the  forward  direction.  The  scalar  theory  of 
Twersky  (1964),  and  Keller  (1964)  for  the  effective  wave  retarded  by  interaction  with  point 
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scattered  does  not  contain  such  a  term. 

SUMMARY 

A  format  method  for  characterization  of  scattered  by  a  moment  expansion  of  their 
equivalent  pseudo  forces  is  defined  in  equations  (10)  and  (11).  The  method  is  shown  to  be 
equivalent  to  the  form  factor,  or  Fourier  transform  approach  in  the  far  field  as  described  in 
equations  (17)  or  (27a, b,c).  The  total  scattered  energy  may  be  obtained  by  integration  of  the 
radiated  strain  energy  (equation  28c),  or  from  the  forward  scattering  theorem  (equation  28b). 
Forward  conversion  of  P-to-S  wave  energy  under  the  Bom  approximation  is  possible  if  the 
scatterer  is  nonsymmetric  or  anisotropic.  The  first  few  terms  of  for  the  Gaussian,  exponential, 
and  ellipsoidal  scattered  are  given  in  equations  (29),  (30),  and  (31)  respectively.  Trade-offs 
between  shape,  inhomogeneity,  and  anisotropy  exist,  but  a  general  scatterer  may  not  be 
modeled  by  a  homogeneous  scatterer  of  arbitrary  shape.  The  randomly  inhomogeneous 
scatterer  is  treated  with  attention  to  the  attenuaton  experienced  by  plane  P  and  S  waves  pro¬ 
pagating  through  a  distribution  of  such  scattered. 
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VI 


CCS  PROGRESS  REPORT 


One  of  the  by-products  of  DARPA  supported  research  at  Berkeley 
over  the  past  few  years  has  been  the  accumulation  of  a  substantial 
library  of  digital  seismograms  from  explosions  and  earthquakes.  These 
data  are  primarily  broadband  and  primarily  recorded  at  near  and 
regional  distances.  Considerable  effort  has  gone  into  archiving  these 
data  so  that  they  might  be  readily  accessible  for  general  discrimination 
research. 

A  major  problem  with  this  seismic  data  library  has  been  that  of 
providing  efficient  computer  access  to  the  data.  A  related  problem  is 
that  of  providing  the  computational  power  necessary  to  analyze  the  data 
once  they  have  been  accessed.  Problems  such  as  these  have  led  us  in  the 
last  year  to  establish  a Computat ional  Center  for  Seismology  (CCS)  at 
Berkeley.  y 

The  concept  of  CCS  began  to  take  shape  in  late  1981  and  developed  to 
the  stage  of  a  joint  proposal  in  June  of  1982.  The  principal  investigators 
were  T.  V.  McEvilly,  S.  Coen,  L.  R.  Johnson,  and  E.  L.  Majer,  who  had  the 
combined  affiliations  of  the  Department  of  Geology  and  Geophysics,  the 
Department  of  Materials  Science  and  Mineral  Engineering,  the  Seismographic 
Station,  and  the  Earth  Sciences  Division  of  Lawrence  Berkeley  Laboratory 
(LBL) .  Start-up  funds  were  provided  by  developmental  funds  of  LBL  and 
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che  Department  of  Energy,  Office  of  Basic  Energy  Science  (DOE-OBES), 
and  che  center  was  operating  by  August,  1982. 

CCS  has  been  established  as  an  organizational  group  within  the 
i'arth  Sciences  Division  of  LBL  and  is  physically  located  in  the  computing 
center  at  LBL.  It  works  closely  with  the  Computer  Science  and  Mathematics 
(CSAM)  group  at  LBL.  CCS  has  been  set  up  so  that  it  has  access  to  all 
of  the  facilities  of  a  large  modern  computer  center  but  does  not  have  to 
deal  with  any  of  the  operational  details  of  running  such  a  computer  center. 
Figure  1  is  a  somewhat  schematic  diagram  showing  the  relationship  between 
CCS  and  the  major  computational  facilities  at  LBL.  Some  of  the  equipment 
specifically  dedicated  to  the  seismological  computational  needs  of  CCS  is 
attached  to  a  particular  VAX  11/780  computer,  but  CCS  personnel  have 
access  to  the  entire  set  of  LBL  computers  shown  in  Figure  1.  In  addition, 
there  is  the  multitude  of  data  storage  devices  and  input/output  devices 
which  are  accessible. 

At  present  about  a  dozen  faculty,  students,  and  LBL  scientists  are 
actively  using  CCS  facilities.  The  DISCO  package  for  processing  seismic 
reflection  data  has  been  installed  as  part  of  the  CCS  software.  A  variety 
of  other  seismological  software  is  in  various  stages  of  being  transported 
to  CCS  or  developed  within  CCS.  Considerable  emphasis  is  being  placed 
on  the  development  of  efficient  software  for  interactive  processing  of 
seismological  data.  Through  interaction  with  CSAM  and  direct  contact 
with  the  Center  of  Seismological  Studies  (CSS)  in  Arlington,  Virginia, 
there  is  a  concerted  effort  to  maintain  a  general  comparability  with 
otner  DARPA  efforts  in  the  area  of  seismological  computing. 


Figure 


